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AN IMPLICIT FINITE -DIFFERENCE SOLUTION TO THE 


VISCOUS SHOCK LAYER, INCLUDING THE EFFECTS 
OF RADIATION AND STRONG BLOWING 

By L. Bernard Garrett, G. Louis Smith, 
and John N. Perkins* 

Langley Research Center 

SUMMARY 

An implicit finite -difference scheme is developed for the fully coupled solution of 
the viscous, radiating stagnation-streamline equations, including strong blowing. Solu- 
tions are presented for both air injection and injection of carbon-phenolic ablation products 
into air at conditions near the peak radiative heating point in an earth entry trajectory 
from interplanetary return missions. A detailed radiative -transport code that accounts 
for the important radiative exchange processes for gaseous mixtures in local thermody- 
namic and chemical equilibrium is utilized in the study. 

With minimum number of assumptions for the initially unknown parameters and 
profile distributions, convergent solutions to the full stagnation -line equations are rapidly 
obtained by a method of successive approximations. Damping of selected profiles is 
required to aid convergence of the solutions for massive blowing. It is shown that certain 
finite -difference approximations to the governing differential equations stabilize and 
improve the solutions. 

Detailed comparisons are made with the numerical results of previous investiga- 
tions. Results of the present study indicate lower radiative heat fluxes at the wall for 
carbon -phenolic ablation than previously predicted. 

INTRODUCTION 

A blunt spacecraft entering planetary atmospheres at earth hyperbolic speeds 
encounters intense radiative heating rates, particularly in the frontal or stagnation region. 
Acceptable interior temperatures can be maintained by mass -transfer cooling through the 
use of heat shields constructed of polymeric ablator materials (Walberg and Sullivan 
(ref. 1)). Near altitudes for peak heating in many entry trajectories, such as the trajec- 
tory for entry into the Earth's atmosphere from. a. direct -return Mars mission, the mass 
of gas from ablation products injected into the flow field is an appreciable fraction of the 
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mass of the oncoming flow (Chin (ref. 2)) and can blow the viscous boundary layer off the 
surface of the spacecraft. 

This condition, which is generally referred to as strong, or massive, blowing, has 
a physically destabilizing effect on the flow field which seriously impairs numerical solu- 
tions to the governing equations (Libby and Sepri (ref. 3)). Libby (ref. 4) describes the 
physics of the problem to be that of an inner region near the wall which is dominated by 
pressure and inertia forces and in which viscous effects are small, and a thin viscous 
region which adjusts to the flow in the outer inviscid region. For small blowing rates, 
the viscous boundary layer is attached to the wall, and the presence of the solid wall has 
a stabilizing effect on the flow. However, for the strong -blowing problem the boundary 
layer is not attached to the wall, and the gaseous inner layer may destabilize the flow. 

The problem of massive blowing where the viscous flow adjusts to the outer inviscid 
flow has been studied extensively for nonradiating flows (cf . , e.g., refs. 3 to 7). However, 
there has been limited attention directed to the solution when radiation is coupled into the 
problem. Numerous previous investigators using a variety of assumptions and numerical 
techniques have arrived at different conclusions regarding the effects of ablation products 
on the attenuation of radiation. These approaches either required excessive computer 
time for the numerical solution or were approximate analyses which led to questionable 
inputs of the thermodynamic properties required for the radiation computations. Since 
radiative heat fluxes (in addition to being strong functions of temperature and density) 
can also be sensitive to the chemical species within the shock layer, some of which are 
strong radiation emitters and absorbers (Hoshizaki and Lasher (ref. 8)), care should be 
taken in defining these quantities across the entire shock layer. 

The purpose of this investigation is to develop an approach for the numerical solu- 
tion to the coupled equations for the viscous, radiating flow field along the stagnation line 
of a blunt body under both weak and strong blowing conditions. The method solves the 
governing Navier-Stokes equations without making unnecessary simplifying assumptions 
to the equations or in the numerical solution which could degrade the radiative -flux 
computations. 

An implicit finite -difference method, which has previously been shown to be com- 
putationally efficient for chemical -nonequilibrium studies without mass addition (Blottner 
(ref. 9)), is developed for the solution to the proposed problem. The nonlinear governing 
differential equations are written in finite -difference form at all nodal points within the 
shock layer, with boundary conditions specified at the wall and immediately behind the 
shock. The formulation results in a tridiagonal matrix system of algebraic equations 
which is efficient for machine computation. The governing equations are solved "one at 
a time” in succession across the entire shock layer. The technique for the overall 
numerical solution to the two -point boundary problem is by iteration of the flow variables 
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at each nodal point in the shock layer. Since the governing equations are solved one at a 
time rather than concurrently at each nodal point (as in initial -value forward -integration 
techniques), then the nodal -spacing requirement for all the equations is not limited by the 
stability requirement of the most unstable equations. 

The stability of the finite -difference solutions to the governing equations is treated 
extensively. Particular attention is directed to the stability of the species -continuity 
equation with binary diffusion by using central differencing and a two-point windward- 
differencing scheme that provides automatic damping of the profiles. 

Radiation computations are carried out with the use of an existing radiative - 
transport computer program, RATRAP, which was developed by Wilson and Hoshizaki 
(e.g., ref. 10) for a nongray gas. The program uses the tangent-slab approximation (one- 
dimensional) that accounts for absorption and emission within a layer of arbitrary optical 
thickness and is for equilibrium gaseous mixtures of hydrogen, carbon, nitrogen, and 
oxygen. Profiles of pressure, enthalpy, and elemental composition are computed by the 
viscous -flow -field solution and are used as inputs to the radiation program. Radiative 
fluxes are evaluated at various nodal points within the shock layer to provide coupling 
with the flow field. The transport properties for equilibrium air (i.e., viscosity and 
reactive Prandtl number, see Hansen (ref. 11)) are used in the computations to account 
for the energy transport due to binary diffusion. 

Numerical results are presented for both air-to-air injection and for the injection 
of the ablation products of a carbon -phenolic ablating heat shield into an airstream for a 
range of blowing rates of interest. Additional computations are presented for constant- 
density flows and for viscous, nonradiating flows with air-to-air injection. 

The air-to-air injection cases include calculations for a body with a spherical nose 
radius of 3.05 meters entering the earth’s atmosphere at a velocity of 15.24 kilometers 
per second, at an altitude of 61 kilometers, and with a ratio of blowing mass rate to free- 
stream mass-flow rate of 0.1. The results are compared with the initial-value forward- 
integration method of Rigdon et al. (ref. 12), the integral approaches of Smith et al. 

(ref. 13), and an "exact" numerical solution of Wilson (ref. 14), for identical free -stream 
and wall boundary conditions. An assessment is made of some approximations contained 
in these analyses. 


REVIEW OF LITERATURE 

Recent reviews of the advances in analyses for the radiating shock layer are giyen 
by Anderson (ref. 15), Garrett (ref. 16), and Goulard et al. (ref. 17). Investigations related 
to the present work are reviewed in this section. 
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Radiation-Transport Models 

Some existing radiation -transport computer programs which consider ablation 
products in addition to air chemical species are RATRAP, developed by Hoshizaki and 
Wilson (see Wilson (ref. 10)); SPECS, developed by Thomas (ref. 18); and RADICAL, 
developed by Nicolet (refs. 19 and 20). These programs, which perform radiation com- 
putations based on the spectral details of the -individual chemical species, are time con- 
suming for stagnation -line analyses. However, usage of this type of program is required 
since realistic analyses have indicated that the addition of ablation products to the flow 
field can reduce the heat fluxes to the body by factors of 2 below those for pure -air flow. 
(See, e.g., Hoshizaki and Lasher (ref. 8), Coleman et al. (ref. 21), and Chin (ref. 2).) 

For air, significant differences still exist not only in the spectral absorption coef- 
ficients for certain chemical species but also in the radiation models that are coded for 
computer solutions (Suttles (ref. 22)). It is noted that differences in the heat -flux pre- 
dictions at the body can result not only from differences in the radiation -transport codes 
but also from the inaccuracies in the solutions to the flow equations. In particular, the 
predicted concentration, or number density, of strongly absorbing or emitting chemical 
species in the ablation products can be greatly affected by computed temperatures and 
densities within the shock layer. 

RATRAP is used in this analysis. This program is somewhat time consuming but 
contains the appropriate detail for the rigorous flow-field computations of the present 
analysis. 


Flow -Field Analyses With Radiation 

The appearance of the divergence of the radiative flux as a term in the general 
energy equation describing the flow of a radiating gas couples the flow -field and the 
radiative -transport analyses. For hyperbolic entry velocities, typical of return missions 
from Mars, the radiative heat fluxes are large, and thus radiative cooling (energy losses) 
with the shock layer must be considered. The net effect of radiative cooling is a reduc- 
tion of the radiative heat flux to the body because of a decrease in temperature and a 
subsequent increase in shock -layer density from that for adiabatic conditions. (See, e.g., 
Wick (ref. 23) and Hoshizaki and Wilson (ref. 24).) This coupling between radiative 
transport and gas dynamics for this problem is commonly referred to as the "coupled 
radiating shock-layer problem." 

Howe and Viegas (ref. 25) were the first to investigate the flow of a viscous, radiat- 
ing, self -absorbing gas in the stagnation region with the effects of mass addition included. 
Since they assumed a gray radiation model, the radiation -flux results are not quantitative. 
For the axisymmetric stagnation region, they reduced the Navier-Stokes equations to 
ordinary differential equations. They solved the momentum equation numerically by an 
initial -value forward -integration method. The energy and species -continuity equations 
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were solved by numerical evaluation of the exponential integrals that appear in the exact 
analytical solutions to the differential equations. The solution to the coupled system of 
equations was iterated by convergence on the value of the shear stress at the wall. In 
the fully viscous analysis of Howe and Viegas the mass injection rates were restricted 
to moderate amounts which did not upset the stability of the boundary layer near the wall. 
The difficulty with the stronger blowing rates can be traced to numerical instability prob- 
lems associated with evaluating the exponential integrals in the exact analytical solutions 
to the governing equation. 

Wilson (ref. 14), who applied an approach similar to that of Howe and Viegas 
(ref. 25) (with the notable exception being his momentum -equation solution), describes 
the problem for the fully viscous shock layer with massive blowing to be one of numeri- 
cal precision required to take the differences between the large numbers which appear in 
the exponential integrals. Further aspects of this problem will be discussed in associa- 
tion with Wilson's paper (ref. 14) and in the section "Results and Discussion." 

Hoshizaki, Wilson, and Lasher developed integral methods for the coupled solution 
of the viscous,. radiating shock layer around a blunt body (refs. 24, 26, and 8). In each 
paper, the velocity profile was expressed as a fifth-order polynomial. In reference 8 by 
Hoshizaki and Lasher, the species continuity and energy equations were solved by means 
of similarity transformations and numerical integration of the resulting exponential inte- 
grals. Their detailed analysis of the spectral absorption coefficients for 20 air and 
ablation -product chemical species showed that the carbon atoms, which diffuse far into 
the shock layer, act as strong radiation absorbers. 

Chin (ref. 2) developed a numerical method for solving the coupled equations for the 
radiating, inviscid stagnation flow with mass addition. He integrated the conservation 
equations for the air layer from the shock wave to the interface between the air layer and 
the ablation -products layer, and from the body to the interface. He iterated on the wall 
heat flux, the blowing rate, and the velocity and enthalpy profiles until the solution con- 
verged on the enthalpy distribution and heat flux to the wall. His solutions to the inviscid 
governing equations converged very rapidly; however, since the viscous region (which 
typically occupies about 10 percent of the shock layer for the Reynolds numbers of inter- 
est) was neglected, there was no mechanism provided for the diffusion of the strongly 
absorbing carbon atoms and ions into the air layer. 

The first numerical solution to the exact. Navier -Stokes equations for the thin shock 
layer at the stagnation line, including radiative transport, was presented by Rigdon, 
Dirling, and Thomas (ref. 27). This analysis included massive blowing (up to 10 percent) 
for air-to-air injection. In 1969 they extended the solution to massive blowing of abla- 
tion products (ref. 12), with blowing rates up to 20 percent. The numerical procedure 
which they employed is computationally time consuming. They used an initial -value 
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forward -integration scheme in which they were required to adjust four initial conditions 
(temperature, tangential velocity gradient, temperature gradient, and the shear) evaluated 
at the stagnation point. These initial conditions were iterated until the two boundary con- 
ditions for both the momentum and energy equations were satisfied at the shock and at the 
body. In the ablation analysis they were further required to satisfy the binary species - 
coiitinuity equation over the shock layer. If the same initial-value forward -integration 
scheme is used, then this requires a guess for the mass fraction of the injected foreign 
material and its gradient at the stagnation point in order to satisfy the boundary condi- 
tions at the shock wave and the body. The net result is the requirement that one guess 
a priori six coupled initial conditions which are unknown in order to satisfy the boundary 
conditions. 

Rigdon et al. (ref. 12), from the results of their solution to the exact governing 
equations, were able to make direct comparisons with the integral results of Wilson and 
Hoshizaki (ref. 28) and concluded that the polynomial approximations which had worked 
so well in the solution for nonblowing were not sufficient to describe the solution for the 
momentum equation in the presence of massive blowing. Although differences existed in 
the radiation models (Rigdon et al. (ref. 12) used the SPECS code of ref. 18, whereas 
Wilson and Hoshizaki (ref. 28) used an updated version of RATRAP (see ref. 10)), the dif- 
ferences by factors of 2 to 4 in the ablation-layer thicknesses could not be explained on 
the basis of differences in radiation models. 

In reference 14 Wilson concluded that the approximate integral solution to the 
momentum equation was inadequate for large mass injection and/or the Reynolds numbers 
of interest. As mentioned previously, his treatment of the energy and species -continuity 
equations was similar to that of Howe and Viegas (ref. 25). He used similarity conditions 
to obtain exact analytic solutions to the energy and species -continuity equations. He 
attributes the momentum -equation solution (see ref. 28) to a solution technique developed 
by Chou and Blake (ref. 29) for a similar problem. Upon performing an additional coor- 
dinate transformation which involves viscosity, and assuming that the density -viscosity 
product is a constant across the shock layer, Wilson (ref. 14) developed the second-order 
differential equation. He subsequently differentiated the equation, which makes it linear 
(in the second derivative), to obtain an exact analytical solution to the equation in terms 
of exponential integrals. Rather than solve the momentum equation with an initial-value 
technique as was done by Howe and Viegas (ref. 25), Wilson employed a successive 
approximation algorithm to all the analytical governing equations. In the successive 
approximation scheme, which is similar to the technique applied in the present study, the 
distribution of properties is initially specified across the shock layer and the governing 
equations are iterated until satisfactory convergence is obtained. 

Wilson (ref. 14) observed that with his formulated equations he was unable to obtain 
a numerical solution to the fully viscous equations for relatively large blowing rates 
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(greater than 5 percent, approximately). As previously mentioned, Wilson traced the 
problem to numerical instability associated with taking the differences between exponen- 
tially large numbers (greater than about e^O) which were about the same order of magni- 
tude. Since the computer carries only about eight to 16 significant digits, then the result- 
ing difference between these large numbers becomes meaningless. Apparently, these 
exponentially large numbers occur in the numerical solution in regions where the viscous 
effects become small (Wilson observed the effect in the inner region near the wall for 
large blowing). Part of the problem could be due to the loss in precision when Simpson's 
rule or the like (Conte (ref. 30)) is used to integrate numerically under an exponential 
curve . 

Wilson (ref. 14) was able to circumvent the numerical precision problem for mas- 
sive blowing by solving inviscid equations in the inner region and the fully viscous equa- 
tions at a somewhat arbitrary distance from the body. The interface criteria was that e 
be raised to a power less than 10. 

In reference 13, Smith, Suttles, Sullivan, and Graves presented a combined flow- 
field and ablation study of a blunt body entering the earth's atmosphere at interplanetary 
return velocities. The analysis yielded transient histories of ablator mass -loss rate for 
a complete entry trajectory. Smith, Suttles, et al. (ref. 13) examined the entire subsonic 
flow region surrounding the blunt body by dividing it into three interacting regions: an 
inviscid outer layer, a boundary layer, and a charring heat shield. The flow in the invis- 
cid outer layer with radiation was determined by a one -strip integral method developed 
by Suttles (ref. 31). The combined solution for the inviscid flow field with ablation pro- 
vided the boundary conditions for the radiating -boundary -layer computations. The 
boundary -layer and the ablation calculations were iterated until the heating rates and the 
ablation rates converged. For the larger blowing rates the numerical method for the 
solution of the boundary -layer equations was not suitable and they adapted an integral 
procedure by Libby (ref. 6) for the radiating -boundary-layer solution. The solutions for 
the boundary layer and the inviscid outer layer were joined by assuming a cubic variation 
for the distribution of the elemental mass fractions within the boundary layer and by 
adjusting the edge enthalpy condition. 

One of the significant conclusions of Smith, Suttles, et al. (ref. 13) was that the 
introduction of ablation products into the boundary layer did little to attenuate the radia- 
tive flux to the wall. In their analysis they used the RATRAP computer code developed 
by Wilson (ref. 10). Whereas Chin's solutions (ref. 2) for a blowing rate of 7.6 percent 
and the solutions of Rigdon et al. (ref. 12) for a blowing rate of 20 percent indicated reduc- 
tions in the wall heat flux of about 45 percent below the nonblowing rates, Smith, Suttles, 
et al. (ref. 13) calculated reductions of only about 22 percent. Wilson and Hoshizaki 
(ref. 28) in their approximate integral approach indicated (by use of the RATRAP radiation 
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code) that for a blowing rate of 10 percent, the radiative heat flux was reduced 40 percent 
below the value for no blowing, and for blowing rate of 20 percent, the heat flux was 
reduced 60 percent below the no -blowing value. Wilson's more recent results (ref. 14), 
using his improved momentum -equation solution, have indicated much lower heat -flux 
reductions (only 18 percent) at blowing rates of 5 and 10 percent. Although there are dif- 
ferences in the free -stream conditions in the studies of references 28 and 14, they do not 
appear to be sufficient to account for the differences in the radiation blockage effects in 
the two analyses. Apparently, the answer to these differences must reside in part in the 
analytical and numerical treatments of the governing flow equations. 

One of the questions which the present analysis will attempt to answer is whether 
the approximate integral treatments of the governing equations and/or "exact" numerical 
treatments of approximate systems of equations can accurately define the flow properties 
required for the radiation computations with and without blowing. The numerical tech- 
nique to be developed for the analysis of the coupled, viscous, radiating flow along the 
stagnation line with massive blowing included is an implicit finite -difference scheme. 
Blottner (ref. 9) has shown this approach to be computationally superior to initial -value 
schemes for studies of the stagnation -line flow with chemical nonequilibrium and without 
blowing. 

In the implicit method, the problem is treated as a two-point boundary -value prob- 
lem in which boundary conditions are specified at the shock and at the body. The entire 
shock layer is treated as viscous; thus, no "patching" of two or more solutions is required. 
The governing equations for the thin shock layer which describe the viscous flow along 
the stagnation line (Ho and Probstein (ref. 32)) and which are exact through second order 
are solved at discrete nodal points along the stagnation line by iteration through the appli- 
cation of a method of successive approximations. Singularities do not appear in the for- 
mulation since the viscous term takes effect as the convective terms (mass flux) approach 
zero in the vicinity of the stagnation point. 

SYMBOLS 


a tangential velocity gradient 


A n coefficient of the (n - l)th unknown quantity for the nth governing equation 

(see, e.g., eqs. (52) and (53)) 

B n coefficient of nth unknown quantity for the nth governing equation (see, e.g., 

eqs. (52) and (53)) 
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Planck function defined by equation (17), watt-sec/cm^-sr 
arbitrary function (see eqs. (42)) 
specific heat 

coefficient of the (n + l)th unknown quantity for the nth governing equation 
(see, e.g., eqs. (52) and (53)) 

damping factor for the enthalpy solution 

density damping factor 

binary diffusion coefficient for the ith and jth chemical species 

known function appearing on the right-hand side of the nth governing equation 
(see, e.g., eqs. (52) and (53)) 

modified known function appearing on the right-hand side of the nth governing 
x-momentum equation (see eq. (60)) 

exponential integral of order 2 (see eq. (17)) 

arbitrary function (see eqs. (42)) 

static enthalpy 

« 

Planck’s constant, # = 6.6256 X 10“*^ j-sec 
total enthalpy 

mass diffusion flux of the ith chemical species 
thermal conductivity 

Boltzmann constant, k* = 1.38054 x 10”^ J/K 
body curvature, K = l^t^ 



Mi 

molecular weight of species 

i, g/gmole 

N 

total number of nodal points across the shock layer 

Npr 

Prandtl number 


N Re 

Reynolds number defined by 

N 

%e- 

N Sc 

Schmidt number 


P 

static pressure 



qp v>qn v>qp v heat * n y-direction due to conduction, diffusion, and radiation, 

' / W “J 

respectively 

r radius measured from axis of symmetry of body (see fig. 1) 

Rb body radius 

S(x) step function, 

S(x) = 0 for x < 0 
S(x) = 1 for x S 0 

T' temperature, K 

u tangential velocity component 

U resultant velocity 

v normal velocity component 

vj diffusion velocity of the ith chemical species 

x distance measured along the body surface (see fig. 1) 

mole fraction of the ith chemical species 

y distance measured normal to the body surface (see fig. 1) 
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mass fraction of the ith chemical species 


<*v 

0 

6 


Atj 

£ :/c* 

«n 

V 

e b 


K 

A 


M 


P 


a ij 


T 


mass fraction of the ith chemical element 
modified linear absorption coefficient, cm - * 
pressure -gradient term (see eq. (14)) 

normalizing parameter in the y to 77 transformation defined by 

- py s 

5 = \ p(y)dy 
J 0 

nodal spacing 

molecular constant for the ith chemical species (see eq. (25)), K 
molecular constant for the ith and jth chemical species (see eq. (25)), K 
error in the ith iteration at the nth nodal point 
dummy variable 

transformed coordinate defined by equation (30) 
local body angle 

scale factor, k = 1 + Ky = 1 + K 6 P 7 ^ 

Jq P 

eigenvalue for the stability analysis 

viscosity (also dummy variable used in exponential term E 2 in eq. (17)) 
density 

collision cross section "of' "species- -i-j -A 

mean collision-cross. section for species i and j, A 
convergence integral 
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a. net volumetric rate of production of species i 

1 

o(M)* reduced collision integral used in equation (25) 
i] 

Subscripts: 

A air products 

F foreign or ablation products 

inj injected species (at 77 * 0) 

n nth nodal point; n = 1 at the body (77 = 0) and n = N at the shock (77 = 1) 

s shock 

w wall 

*> free -stream conditions 

Superscript: 
i ith iteration 


Primed quantities are dimensional. 

ANALYSIS 


Thin -Shock-Layer Equations 

The governing equations for the steady-state flow of a viscous radiating gas in the 
stagnation region of a blunt axisymmetric body at moderate -to -high Reynolds numbers 
are given by Ho and Probstein (ref. 32) and by Scala (ref. 33). These equations are 

Continuity: 

-^r(r'p’u’) + A.(K’r’pV) = 0 (1) 

x -momentum: 


pV — + n'p'v' — + K'pVv’ = - -^1 + — (kV $2? 
ax' ay' H ax' ay' \ ay', 


( 2 ) 
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y -momentum: 


p. u . + K'pV ^ - K’ p’(u’) 2 = -k’ -|1. 


Energy: 


— " K '' >,v ’ r ’ - apf’ 1 " («C,y + «D,y + 4 R,y) 


r P u aT + 


8 ( , , , , 8u f 

+ T7 k r |l U' -r-r 

ay* V K ay* 


Species continuity: 


(3) 


(4) 


-^-(r'p'u’Q'i) + ^(K'r'pVai) = - 

(i = 1, . . M species) (5) 

Primed symbols are used to denote dimensional quantities; unprimed, dimensionless 
quantities. 

For the body -oriented coordinate system shown in figure 1, the quantity k* , which 
is the coordinate stretching function, is defined by 


k* = 1 + K'y’ 

(6a) 

where 


K* = K ’(x) = i 

R b 

(6b) 

and r satisfies the equation 


dr’ = k’ sin 0^ dx* + cos 0^ dy’ 

(6c) 


In equations (4) and (5), the quantities q£, , q^ y, and q^ are the heat fluxes 

in the y* -direction due to conduction, diffusion, and radiation, respectively; jJ y is the 
mass diffusion flux in the y* -direction; and is the net rate of production of the ith 

chemical species. 

In -comparison-with.the.heat and jnass diffusion fluxes in the y* -direction, the cor- 
responding fluxes in the x* -direction are generally considered negligible in the~stagnation 

region. (See Ho and Probstein (ref. 32).) These fluxes in the x' -direction are assumed 

to be negligible in the present analysis also. -The equations for the thin shock layer are - 

the simplified boundary -layer equations (Navier -Stokes) including the curvature terms 
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for the stagnation region. The equations are considered to be accurate to the order 


(ref. 32), where Npj e = °° ,°° — , when radiation is not considered. 


M 


Stagnation-Streamline Equations 

At the stagnation line (x = 0), the conservation equations for the thin shock layer 
reduce to ordinary differential equations upon expanding the flow variables in the follow- 
ing power series (based on symmetry considerations): 


u = a(y)x + a 1 (y)x d + . 
v = v(y) + Vj(y)x 2 + . . 

P = p(y) + P^y)* 2 + • • 

H = H(y) + H 1 (y)x 2 + . 

«i = <*i(y) + « n (y)x 2 + 




(7) 


J 


applying the geometry relations (eqs. (6)), and then taking the limit as x — 0. The limit- 
ing forms of the governing equations for the stagnation -line flow become 


Continuity: 

d(/c 2 pv) _ 

dy 

x-momentum: 


2/cpa 


( 8 ) 


+ KPV^ + pa 2 + Kpva = /3 
N Re d y V dy/ dy 


(9) 


y -momentum: 


pv^= 

P dy dy 


Energy: 

_ K2 («C,y * «D,y + «R,y) 

Species continuity: 


dy dy 


2 da i 

KTfSV — 

dy 



K 2 U)j 


( 10 ) 


(ID 


( 12 ) 
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where 




where the subscripts °° and s refer to free -stream and shock conditions, respectively. 

Restrictions and Assumptions 

The governing stagnation-line equations are general equations and are restricted 
only by the requirements that P s ^Nr 6 « 1 and the net radiative heat fluxes in the 
x-direction are comparatively negligible. The following basic restrictions and assump- 
tions are imposed on the governing equations: 

(1) The^as is'in local thermodynamic and chem ical equilibrium. 

(2) Diffusion of the chemical species is governed by a binary diffusion model. 

(3) Radiative energy transport occurs within a one -dimensional, infinite, planar 
slab (tangent slab). 
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In order to uncouple the stagnation-line solution from the remaining subsonic field, 
it is necessary to assume the relationship between the shock and the body curvatures at 
x = 0. In the results of this analysis, the shock and body are assumed concentric unless 
otherwise specified. The results of Suttles* analysis for an inviscid, radiating flow field 
(ref. 31) indicate that the assumption is reasonable. 

Preheating of the ambient air upstream of the shock wave due to radiation transport 
(precursor effects) is neglected. Smith (ref. 34), Hoshizaki and Lasher (ref. 8), and 
Rigdon et al. (refs. 12 and 27) indicate that precursor effects begin to become significant 
at velocities around 17 to 18 km/sec and above. Thus, the present solutions, which 
employ the Rankine -Hugoniot conditions (see Hayes and Probstein (ref. 35)) for the dis- 
continuous jump conditions across the shock, will be restricted to velocities somewhat 
lower than this. 


In the numerical solutions, unless otherwise specified, a Newtonian pressure dis- 
tribution (ref. 35) is used to evaluate /3, the term in the x -momentum equation; 

\9x 2 /x=0 

that is, (-^H = -2.0. It should be noted that Wilson (ref. 14) used a value of -3.0, 

W 2 /x=0 

which leads to a thinner shock layer in his calculations. 


The radiation-transport computer code which is used in the radiation computations 
is RATRAP, developed by Hoshizaki and Wilson (refs. 24 and 28 and appendix A of this 
paper). RATRAP considers most of the primary radiating chemical species associated 
with carbon, oxygen, nitrogen, and hydrogen mixtures. The calculations of the detailed 
thermodynamic and chemical compositions for chemical equilibrium are performed in 
the computer code FEMP, developed by Browne et al. (ref. 36). FEMP is included by 
Wilson as an integral part of RATRAP . 


In the analysis, it is assumed that the transport properties for air developed by 
Hansen (ref. 11) apply for both air-to-air injection and injection of ablation products into 
air. This assumption for the ablation products will be superficially analyzed by pertur- 
bations in the pertinent transport property (Prandtl number) to determine its effect. It 
should be noted that Rigdon et al. (ref. 12) ran two identical cases with the exception of 
pure -air transport properties in one and the combined transport properties of ablation 
products and air in the other, and the resulting differences in the radiative heat fluxes 
to the wall were less than 1 percent. 


The continuity and y -momentum equations (eqs. (8) and (10)) are unchanged when 
the foregoing assumptions are applied. The only perturbation in the x-momentum equa- 


tion (eq. (9)) is the 


vax 2 /x=i 


, or j3, term arising from the Newtonian (or specified) 
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pressure distribution in which (—^-\ will be taken to be a constant in y in this 

W 2 /x=0 

analysis. 

By applying the binary -diffusion and tangent -slab approximations, the heat -flux 
terms in the energy equation (eq. (11)) become: 


q C,y + q b,y 



(16) 


where k' is the "lumped” or "effective" conductivity computed and tabulated by Hansen 
(ref. 11) for air at temperatures up to 15 000 K and pressures up to 100 atmospheres 
(1 atm = 101 325 N/m2). 

The radiative -heat -flux equation for the heat flux at a point within a one -dimensional 
slab is given by Suttles (refs. 31 and 37) and is 



(17) 



is the exponential defined by 



dji 


and is described by Kourganoff (ref. 38). In equation (17), B u is the Planck function, 
v is the frequency of radiation, and a v is the modified linear absorption coefficient 
which is a function of temperature, pressure, and the number densities of the chemical 
species within the slab. The equation is valid for a nongray self -absorbing gas. The 
absorption coefficients and radiative heat fluxes are computed in the radiation computer 
code RATRAP developed by Hoshizaki and Wilson (see Wilson (ref. 10)). The radiation 
model takes into account both continuum and atomic-line radiation exchange in a slab of 
nonuniform temperature. Details of the radiation model are given in appendix A. 

The simplified energy equation can now be written as 


2 , , dH f 
K p V df 



(18) 
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where q™ v is given by equation (17). The combined conduction-diffusion heat -flux 

^ >y 

term can be written in terms of total enthalpy from the following relations for static 
enthalpy, total enthalpy, and "effective" Prandtl number: 


dh' = Cp dT’ 



( 19 ) 


N Pr = 



( 20 ) 


to yield 


* 2 pV dH 1 = JL| 


kV /dH' _ v , dv’ 
Npr \dy' dy* 




dy’ dy' 

Upon applying the nondimensional relations defined in equations (15) and rearranging 
equation (21), the dimensionless energy equation becomes 


( 21 ) 


k 2 dH = _1_ d/K^v^l + d_/2 \ 

N Re dy \Np r dy / dy N Re dy\N Pr dyy dyV R ,YJ 


( 22 ) 


The species -continuity equation (12) undergoes considerable modification because 
of the assumptions of binary diffusion and equilibrium chemistry. The species -continuity 
equation was given by equation (12) as 


/c 2 pv 

j. 


T5T--&( A ‘jr) 


dy 


K 2 d>i 


Under the assumption of local chemical equilibrium, the volumetric rate of produc- 
tion of species i is indeterminate. However, the fact that the net rate of production of 
the mass fractions of the chemical elements is zero can be utilized to yield a tractable 
solution to the species -continuity equation. If it is assumed that the binary -diffusion law 
holds for the mass fractions, that is, 


_ dQ!j 

J i - -<> d 12 of 


(23) 


where the barred quantities refer to elemental mass fractions, the species -continuity and 
diffusion equations can be combined to yield the elemental -diffusion equation 




dofj 




dy 


dy 


dy 


(24) 


The dimensional binary diffusion coefficient is given by Hirschfelder et al. (ref. 39) 


as 
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r ~ » -.1/2''' 
(T* )' i (M i + Mj) 


Dl. = 0.002628 


L 


2M i Mj 


J 




p’crl.fi^’ 1 ) 
1J. 1J 


J 


where 


a.’. = - (o'. + a!\ 

i] 2V 1 ]J 


(25) 


and 





1/2 


The molecular constants cr^ and e^/k* are tabulated by Svehla (ref. 40) and are given 

(1 11 * 

in table 1 for the chemical species of interest. The reduced collision integral ’ 


is based on the Lennard -Jones potential and is defined in reference 39 as a function of the 
nondimensional term T? 1 . which is given by 



For the binary diffusion model, the individual elements of the ablation products and 
of the air which passes through the shock are considered to diffuse in the same respec- 
tive manner as the two chemical species (one for the ablation products, the other for the 
air) used to form the binary diffusion coefficient. Thus, a distinction need not be made 
between the individual elements but only between the total mass fractions which represent 
the ablation products and the remaining mass fraction which represents the air products. 
Since the total mass fraction of all the elements (and, for that matter, the chemical 
species) must equal unity at any point, then one need solve only one elemental -diffusion 
equation, which is given by 


d /„2 r> dQ! F N ' 

^7 r pD 12 


dy 


dy 


- K^ftV 


-da F ~ 


dy 


= 0 


where ap is the total mass fraction of.the elements of the ablation products, 
mass fraction of the air products a ^ is then computed by 



(26) 
The total 


(27) 
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( 28 ) 


The mass fractions of the individual elements are then calculated by the equations 


a i “F^ipnj 


+ at .a. 

A 1,00 


(i = 1, . . M) 


where ^ is the specified elemental mass fraction of the ith element of the ablation 
products injected at the wall and 3^ is the same ith element that passes through the 
shock layer. The density and the individual chemical species are then calculated from 
an equation of state by the free -energy -minimization subroutine FEMP (ref. 36) where 



(29) 


For analyses of the viscous, radiating stagnation -line flow, the resulting governing 
equations ( 8 ), (9), (10), (22), and (24) are the most general system of equations that are 
treated in the literature. These are the exact equations for the thin shock layer (for the 
tangent-slab and binary -diffusion assumptions) which are solved in this analysis for the 
stagnation line. As previously mentioned, Rigdon et al. (refs. 12 and 27) have applied 
initial -value techniques for the numerical solution to this system of equations (with the 
exceptions that they assume constant pressure across the shock layer and neglect the 
curvature terms) for the massive -blowing problem. 


Transformed Equations 

Before the solution technique is developed, it is desirable to transform the govern- 
ing equations by a change in independent variable from y to 17 , where 

v = l P p(y)dy (30) 

6 J 0 

where 6 is the normalizing parameter given by 

6 = f Ys p(y)dy 

J 0 

This transformation has the important effect of fixing the shock boundary at 17 = 1. The 
inclusion of density in the transformation gives, for a fixed nodal spacing in rj, a finer 
resolution on a physical scale y at points near the body and within the boundary layer, 
where conditions are rapidly changing. 

The transformed system of equations become: 


Continuity: 



x-momentum: 


_se_if 4P fe 

N Re 6 2dT 'V *i 


— - pa^ + Kpva = 
6 dTj 


-/3 


y -momentum: 


^P = -pv^I 
drj d7] 


(32) 


(33) 


Energy: 


d / fc 2 pp dH^ 

d 7 ? \Np r d77y 


£ XT 2 dH d / k^u dv\ j XT d / 

+ 6N Re K p v — = — (rr-* 1 p v — I + 6n rp — (« <1 


dtj d? 7 \Np r drj 


Re 


d vV H *,v 


Elemental diffusion: 


d [ 2 2 n 

Tn (* p D i2 


daj 

drj 


+ dK^pV 


doj 

d?7 


= 0 


(34) 


(35) 


Equation of state: 

p = p(a v h,p) (36) 

Note that the normal velocity v has been redefined as negative in the positive y- 
or ??- direction in equations (31) to (36) . This change in no way affects the solution to the 
equations since the boundary conditions, as will be discussed subsequently, have been 
appropriately modified to reflect this change. The flow -field coordinate system for the 
stagnation line in the transformed coordinates is shown in figure 2. 

The transformed coordinate r) is specified as zero at the body and as unity at the 
shock. The normalizing parameter 6 is an unknown in the problem and is obtained 
from the solution of the continuity equation (31): 

s = v<4» - (f 2p ik m 

2 \ /ca drj 

J 0 


Boundary Conditions 


The subsonic flow field in the nose region is described by elliptic equations; conse- 
quently, the stagnation-streamline'solutTdrfis’ihfluence'd'by'th'e'flow within-the-entire-sub- 


sonic region. However, this influence only enters in the 0 term, that is, 


.m 


,ax 2 / x =o 


in the x-momentum equation (32), and in the curvature of the shock wave which provides 
the boundary condition for the tangential velocity gradient at the shock wave. For a shock 
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wave and body which are concentric, the Rankine-Hugoniot relations 

/ q2J\ 

and the velocity gradient a s (see ref. 


relation between 


. 9x2/x=0 


yield the following 
35): 


0 = - 





Newtonian impact theory predicts (ref. 35) 


v 9 x 2 /x=0 


- 2.0 


Thus, 



1.0 


(38) 


(39) 


With /3 and a s specified, the stagnation-streamline problem becomes a two- 
point boundary -value problem, where conditions are specified at the shock wave and at 
the body. The boundary conditions necessary to solve the transformed conservation equa- 
tions for the flow of a viscous, radiating gas with injection of foreign species, where the 
chemical species are in local thermodynamic equilibrium, are 

At the body surface, or wall (77 = 0, y = 0): 


a = 

II 

O 


pv 

II 

? 


P = 

Pw“P S 

+ |(p v2 )s 

H = 

= H W 




(and a-p 




> (40) 


(i = 1 , . . ., M elements)J 


where the boundary condition Op = 1 applies for conditions of strong blowing. 
At the shock wave (77 = 1, y = y s ): 
a = a s = 1.0 
pv = (pv) s = 1.0 

P = P S \ 

H = H s = H to 

cfj = ®i,oo ( an d Op = 0 j (i = 1, . . ., M elements) 


(41) 
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where Q' i ini is the mass fraction of element i which is injected at the wall, and 
0 !^ is the mass fraction of element i which passes through the shock wave from the 
free stream. Conditions immediately behind the normal shock are computed from 
Rankine-Hugoniot relations. Heat fluxes entering the slab boundaries are assumed to be 
zero; that is, there is no precursor heating of the free stream by the shock layer and no 
radiation from the body into the shock layer. 


Numerical Solution 

The finite -difference approximations to the governing equations (31) to (36) and the 
procedure for the numerical solution to this system of equations are developed in this 
section. The equations are written in finite -difference form for a network of N equally 
spaced (in rj) nodal points between the body (n = 1) and the shock (n = N) , which are shown 
in figure 3. For most of the calculations, 21 nodal points are used. The overall numeri- 
cal solution technique is iteration. Properties and parameters are assumed initially for 
each nodal point across the shock layer. The governing finite -difference equations are 
then solved sequentially, with the use of the most recent values of the profile parameters, 
until satisfactory convergence is obtained at each nodal point. 

Finite -difference approximations . - Either two- or three-point finite -difference for- 
mulas are used to differentiate and integrate the governing differential equations numer- 
ically. Whenever possible, three -point central differences are employed since they 

p 

are accurate (at a particular nodal point) to order (At?) , where A 77 is the distance 
between points; whereas two -point differences are accurate only to order A77. (See, e.g., 
ref. 30 or 41.) Two-point windward differences (with the flow) are employed for the con- 
vective terms in the elemental -diffusion and energy equations where dictated by stability 
requirements (i.e., profile solutions which are well behaved across the shock layer). 

The stability study of these two equations is presented in appendix B. 

The finite -difference formulas that are employed to approximate the derivatives 
are given below. (See, e.g., ref. 30 or 41 for the development.) The formulas are valid 
for equally spaced increments in At?. The central-difference formulas are 


I df\ _ ~^n-l + ^n+1 
\dr?j n 2 A 77 


(42a) 


f d 2 f\ _ *n-l 2 ^n + *n+l 
V d7 ? 2 /n (At ?) 2 

d _ _c n f n-l " ( c n-l ~ c n+l) f n + c n f n +l 

_dT? ( } J n = 2 At? 


(42b) 


(42c) 
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d / df \ _ ( c n-l + c n)*n-l ~ ( c n-l + 2c n + c n+l) f n + ( c n + c n+l) f n+l 

dn [ C d ifj] n " 2 (At?) 2 


(42d) 


where c and f are arbitrary functions. In the development of equation (42d), it was 
assumed that c varies linearly between nodal points. 

The windward -difference formulas are either forward or backward differences, 
depending on the direction of the flow. If the mass flux pv is positive (toward the body), 
forward differences are taken for the convective terms; if pv is negative, backward dif- 
ferences are taken. The formulas are 


If (pv) n > 0, 


pvS.) =(pv) n 


*V n 


If (pv) n < 0, 


f n+l ~ f n 
A77 


*n " *n-l 


pV |) n =(PV)n A 77 


(43) 


(44) 


Numerical integration is performed by the Simpson’s rule approximation (ref. 41) 
using three points at a time along the profiles. The finite -difference equations used in 
performing the quadrature over the interval from Tj n to V n+ \ and over the full interval 
from ?? n to ? 7 n+ 2 are 


and 






(45) 


j" + = + 4I n+1 + f n+2 ) 

% J 

/ 

Initial profiles.- For the general case, the following profiles are initially assumed: 


pv( 77 ) = (pv) w 

+ [(pv) s - (pv) w J? 7 2 

(46a) 

pfa) = Pg + §( 

pv 2 ) s (l - 77 2 ) 

(46b) 


a (t?) = a s 77 


(46c) 


h(v) = h w + (h s - h w )7? 


(46d) 
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(46 e) 


a pfo> = (« F ) W U - V) 


p(v) = pf Q! i( T ?) ,pO?),h(tj)] 


( 461 ) 


The density is computed in subroutine FEMP for a gas mixture in chemical equilibrium. 
All other quantities appearing on the right-hand side of the initial -profile equations are 
either available from the input boundary conditions or they are computed from Rankine- 
Hugoniot relations for the normal shock. 


An initial value for the transformed shock displacement 6 is also required. The 
physical shock displacement distance for no blowing is correlated by Inouye (ref. 42) to be 


yk 


P’oc 


- 4 - * 0.78 -4 


R 


b 


Ps 


Since = Sp' (* then for a constant shock-layer density of 
R b °° J 0 P 


'0 P 
6 * 0.78 


For moderate blowing, 

6 * 1 (47) 

which is the initial value assumed in this analysis. 

This completes the statement of the problem of determining the initial profiles for 
the general case. The appropriate boundary conditions and initial profiles for special 
cases of interest, such as no blowing or air-to-air injection, are covered in the subse- 
quent discussion of these special cases. 

Solution procedure and finite -difference equations .- With the assumed initial pro- 
files and parameters (eqs. (46) and (47)), the conservation equations are solved by suc- 
cessive iteration. The coupled equations are solved in the following sequence: continuity, 
y -momentum, elemental continuity, x -momentum, energy, and equation of state. The 
most recently computed values of the profile parameters and 6 are used in the compu- 
tations. The density distribution is necessary to solve the governing differential equa- 
tions. After these equations have been solved, updated values of density are computed 
from the equation of state and compared with the previous density values at each nodal 
point. The entire sequence through the governing equations is repeated with the new 
density distribution until two successive passes yield nearly identical (within 1 or 2 per- 
cent) density values at corresponding nodal points. The flow diagram which illustrates 
the solution procedure is shown in figure 4. 
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With the distribution for the tangential velocity gradient and the scale factor k(tj), 
the continuity equation (31) is numerically integrated by Simpson's rule (eq. (45)) to give 
the updated shock displacement distance 6: 


A 

8 = 


~(pv) w + (*V) S 

2 f ksl dri 
J 0 


(48) 


and the mass -flux distribution pv(ij): 



*a d?? + (pv) w 


( 49 ) 


The velocity profile, which is required in the solution of the y-momentum equation, 
is computed from the previous density distribution p(n) and the updated mass -flux dis- 
tribution pv(??), namely, 


v<„) = 

p(v) 


(50) 


The y-momentum equation (33) is numerically integrated by Simpson's rule (eq. (45)) 
to give the pressure distribution p{rj): 



pv dn + p 

dr/ b 


(51) 


The profiles of pv(r]) and v(?j) are tabulated from the continuity -equation solu- 
tion, and a central-difference numerical differentiation scheme (eq. (42a)) gives 

d M dl1 

at the N - 2 points within the shock layer. The values of =21111 at the body and the 

dri 

shock boundaries are obtained from two-point forward- and backward-difference approxi- 
mations, respectively. 


The elemental-diffusion equation (35) is cast in finite -difference form by using the 
central-difference formula (eq. (42d)) for the second -derivative diffusion term and the 
windward -difference equations (eqs. (43) and (44)) for the first -derivative convective term. 
The resulting equations are 


For n = 1, 

Mi = 1 


(52a) 
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For n = 2, . . N - 1, 

(/c 2 p 2 D 12 ) n i + ( K ^P^12) n " 26(Ar])(/< 2 pv) n (l - S) ( a y) n=1 _ ( K p D 12)p_i 
+ 2( K 2 p 2 D 12 ) n + (« 2 P 2 D l2 ) n+1 + 26(a.l)( K 2 pv)„(l - 2S)](a F ) n 
+ [(« 2 p 2 d 12 )„ + (« 2 P 2 Dl 2 ) n+ i + 2«<4n)(K 2 Pv)„s](3 F ) n+1 - 0 < 5! 


For n = N, 

( 3 f)n = 0 


where 

S = 0 if (pv) n < 0 


(52c) 


S = 1 if (pv) n S 0 

The formulation results in a banded (tridiagonal) matrix equation of the form 
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where A n , B n , and C n are the coefficients of the (^ <y --p) n j» ( a r) n ’ and (“f) j 

terms, respectively, and D n are the terms appearing on the right-hand side in equa- 
tions (52). At the wall boundary n = 1, the coefficients are Bj = 1, Cj = 0, and Dj = 1. 
At the shock boundary n = N, the coefficients are A^ = 0, B^ = 1, and = 0. The 
system of equations is easily solved by Potters' method, which is a form of Gaussian 
elimination that is efficient for solving a banded matrix system of equations. (See refs. 43 
and 30.) The solution to equations (52) yields the distribution of the total mass fraction 
of the ablation products at N nodal points subject to the constraint that 0 ^ g 1. 

This constraint was imposed only during the intermediate iterations and was not neces- 
sary for the final solution. The mass fractions of the individual air and ablation -product 
elements are then computed from equations (27) and (28). 


The transport properties /i and Np r required in the solutions to the x-momentum 
and energy equations are obtained from a tabular lookup as functions of the temperature, 
calculated in the FEMP subroutine, and the pressure, obtained in the solution of the 
y -momentum equation. 

The x-momentum equation (32) is next solved. The equation is cast in finite - 
difference form by using the central-difference expressions given in equations (42) to 
approximate the first- and second-order derivatives. Since the equation is nonlinear in 
the velocity gradient a, an iterative procedure is required for its solution. A quasi - 
linearization approach (see ref. 44) is employed to generate the solution by a rapidly con- 
verging iteration. 


The x-momentum equation is quasi -linearized by the following technique: At a 
given nodal point, let 





(54) 


where the superscripts i and i-1 refer to the values of a ( t?) at the iterations i and 
i-1, respectively. Upon assuming that 

limja^ - a^ 1-1 ^] 2 - 0 (55) 

i— oo 

and x-momentum equation becomes linear in a^(77), where 


^W^aM-ta * 1 - 1 )] 2 as i 


(56) 


r 
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Combining the finite -difference formulas for the derivatives (eqs. (42)) with equation (56) 
yields the finite -difference form of the x -momentum equation: 

For n = 1, 

= 0 (57a) 

For n = 2, . . N - 1, 



26(At?) 



- (Kpv) n + 


(KP) n (W) n _i + 2(pp) n + (pp) 
23 2 N Re (AT]) 2 




r (*P 2v )n (Kp) n [(pp) n + (PM)n+l] 
26(At?) 




25"N Re (Ar])^ 



= /3 + 



(57b) 


For n = N, 


a (0 _ a 

a N " a s 


(57c) 


The subscript n refers to the nth nodal point in the shock layer and the superscript i 
refers to the ith iterative value of a. 


For a linear iteration on the velocity gradient a, let 

= a?- 1 * + 


(58) 


~n ~n ' ~n 

where e n is the error in the ith iteration at the nth nodal point. Substituting equation (58) 


for a^ into equation (57b) yields a system of equations of the form 


■^n e n-l + ®n e n + ^n e n+l 


(n = 2, . . N - 1) 


(59) 


,<i>. -to 


,(0 


where A n , B n , and Ch~are the-coefficients of-the terms- a^j.,_a^ ; ,_and_a^p 
respectively, of equation (57b) for n = 2, . . N - 1. The term D n is given by 


D n = |3 + p 


,(i~l) 


'n 


A a^" 1 ) R a^ i_ 0 _ r a ( i_1 ) 
A n a n-1 B n a n L n a n+1 


(60) 
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for n = 2, . . ., N - 1. Since the boundary conditions are specified at the body and at 
the shock, then the equations at the boundaries become 

(or = Ci = Dj = 0 and Bj = ljj 

> (61) 

'or A n = Cn = D n = 0 and Bj<f = 1) I 

The system of equations (59) and (61) may be written as a tridiagonal matrix equation, 
which is solved by Potters’ method. The velocity gradient a n is updated by equa- 
tion (58) and the process repeated until 

| e n | = t (all n) 

In this study, this assigned convergence interval t is 0.01. 

The total -enthalpy distribution is obtained from the tridiagonal -matrix solution to 
the energy equation (34) in a single pass. As previously discussed, the total -enthalpy 
derivative for the combined convective and diffusive heat -flux terms is included explicitly 
in the matrix solution for the H n terms, whereas the radiative -flux term is computed 
by using the enthalpy profile of the previous iteration and the updated profiles of pressure 
and elemental mass fraction from the solutions for the y-momentum and species -diffusion 
equations, respectively. Consistent with the treatment of the elemental -diffusion equation, 
the windward-difference equations (eqs. (43) and (44)) are used to approximate the convec- 
tive term in the energy equation. Three-point central -difference formulas are applied to 
the remaining derivatives to yield the following finite -difference form of the energy 
equation: 

For n = 1, 

Hj = H w (62a) 

For n = 2, . . ., N - 1, 



2 (At?)' 


+ ( «?££ 

^ N Pr/ n _i \Np r 
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For n = N, 


(62c) 


% = H oo 

where 

S = 0 if (pv) n < 0 

and 

S = 1 if (pv) n 5 o 


The derivative terms appearing on the right-hand side of equation (62b) are evaluated by 
the central-differences formulas (eqs. (42)). 

The density distribution which was required in the solution to the governing differ- 
ential equations can now be updated by the equation of state 

p n = ^[(“On^n^nj (63) 


at each nodal point. The updated values of p n provide the mechanism for iterating the 
governing equations. Upon setting p n equal to either the newly computed values pj^ 

or some fraction of the old values of density plus some fraction of the new values 

( in order to speed convergence of the flow -field solution), the entire procedure, 
beginning with the continuity equation, can be repeated until 





^ e 


(all n) 


(64) 


The interval of convergence e used in this analysis is 0.02 or less. 

Salient features of the implicit finite -difference algorithm .- The implicit -finite dif- 
ference algorithm developed in the previous section is quite simple, yet sufficiently flex- 
ible to treat the problem of viscous, radiating shock layers. 

With the two- and three -point difference approximations for equally spaced incre- 
ments, the governing equations for the thin shock layer can be made amenable to numeri- 
cal solution without making unduly restrictive assumptions for the purpose of yielding 
analytically-tractable .solutions. The tridiagonal -matrix system of equations which results 
can be efficiently solved by Potters’ method7which requires only about—3n._ computations, 
as contrasted with the (n) 2 computations required for inversion of a full matrix. 

Simple linear and quadratic profiles 'can be initially assumed as functions of the 
boundary conditions. Or, one may take advantage of prior knowledge of the solution 
behavior and begin with improved estimates of the profile parameters (i.e., read in the 
parameter values at each nodal point). 
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No singularities appear in the governing finite -difference equations since division 
by pv is not required in the present formulation of the problem. In the region where 
pv (or v) approaches zero, the diffusion terms become important in the elemental- 
diffusion equation (eqs. (52)). Likewise, the viscous terms predominate in;the 
x-momentum and energy equations (eqs. (57) and (62)) near pv = 0. 

The method is an implicit scheme in which the unknown quantities at point n are 
calculated as functions of conditions at surrounding points as well as conditions at the 
point itself. This is in contrast with explicit schemes in which the unknown quantities at 
a point are evaluated solely as a function of conditions at a former point (such as the func- 
tion and its derivative being evaluated at the (n-l)th point and then extrapolated to the nth 
point to determine the unknown function). For a given step size, implicit schemes are 
generally stable (bounded), whereas explicit schemes may or may not be stable. 

In the present approach to the two-point boundary -value problem, the boundary con- 
ditions are specified at the shock wave and at the body. Since the matrix system of 
equations (for a particular governing equation) is completely coupled across the entire 
shock layer (e.g., note the appearance of i- n the (n-l)th, nth, and (n+l)th elemental- 

diffusion equations (eqs. (52)) and the known boundary conditions are included in the sys- 
tem of equations, then the computed distributions between each endpoint are bounded and 
are monotonic in behavior (with the exception of the pressure, which reaches a maximum 
at the stagnation point). Thus, the analyst has near -maximum information about the 
behavior and levels of the unknowns which are to be computed as well as a high degree of 
assurance that the computer run will not be aborted because of overflow (numbers larger 
than the computer will accept). 

Perhaps one of the biggest advantages of the present approach over initial -value 
methods resides in the fact that only a limited number of unknown quantities need to be 
evaluated in order, first, to satisfy the governing differential equations and, second, to 
provide the necessary inputs for the radiative -flux computations. In the present method, 
the only unknowns are the properties themselves (pv, p, etc.) at each nodal point about 
which one has maximum information as to the bounds on the values and the general 
behavior of the properties across the shock layer. In contrast, consider the computations 
which must be performed in an initial -value treatment of the problem (such as Runge- 
Kutta forward integration). In the initial -value approach for the solution of a single 
governing equation, the unknown property and also its derivative must be computed. 
Generally, little information is available as to the behavior of the derivative across the 
shock layer. As a consequence, the function or derivative changes from point to point 
are generally closely controlled, by restricting the step size, to maintain stability as the 
calculations proceed downstream. 
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RESULTS AND DISCUSSION 


In this section the numerical results obtained from the solution to the governing 
equations developed in the previous section are presented and discussed. The solutions 
which were obtained in the various phases in the development of the method are presented 
in the sequence of their development, beginning with the constant -density solutions and 
continuing through the solutions for the viscous, radiating shock layer with injection of 
ablation products included. Results are compared with results from existing approaches 
to the solution of the equations for the radiating flow field. 

All numerical solutions were generated on a Control Data Corporation 6600 digital 
computer. 


Constant -Density Solutions 

The results of the constant -density study are presented in this section. The con- 
vergence behavior of the solution for the reduced system of equations that describe the 
flow of a constant -density gas is examined. Selected solutions are shown for inviscid 
and viscous flows. Results are presented for a range of blowing rates and Reynolds 
numbers. 


The overall solution to the governing equations is by iteration until satisfactory con- 
vergence is obtained on the density. For a constant -density assumption (and a viscosity 
and Reynolds number specification required for the viscous -flow solution), the continuity, 
x-momentum, and y-momentum equations can be solved independently of the energy and 
diffusion equations. Naturally, the density is not iterated; however, the remaining sys- 
tem of equations (continuity and x- and y-momentum) are still coupled, and in the present 
approach, an iterative scheme is required for their solution, as indicated by figure 4. 

The equations are solved "one at a time" by a method of successive approximations until 
the computed values of pv, a, and p for two consecutive iterations are within a speci- 
fied accuracy at each nodal point. The constant -density solutions thus serve a twofold 
purpose. First, the solutions provide a fundamental means of studying the convergence 
behavior of the coupled system of equations for a wide variety of blowing rates and 
Reynolds numbers and thus permit an assessment of the adequacy of the successive - 
approximation iterative scheme. Second, the solutions provide an indication as to the 
nodal-spacing r^equirefnents for-the formulation of the implicit finite -difference equations. 


Several cases were run for. a range of typical blowing rates ((pvj w = 0 tcT -0:2\; — - 
Reynolds numbers ^Nj^ e = 1 to loA), and nodal spacings (N = 11 to 101) for a constant 
density of 20 and a viscosity of 1.0. All solutions required about three-or .four iterations 
and a total computer time of 2 to 4 seconds. No discernible differences were observed 
in the computed profiles of pv, p, and a or the values of 6 for N = 11 and N = 101. 
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It was observed that the computer times were a function of the number of nodal points 
(the 2- and 4 -second run times were for N = 11 and N = 101, respectively), and increas- 
ing the blowing rate or Reynolds number had little or no effect on the computer run times. 


The convergence behavior of a typical constant -density solution is shown in figure 5, 
where the successive solutions to the continuity and x-momentum equations are plotted. 
The solution is for a constant density of 20, a constant viscosity of 1.0, a typical Reynolds 
number of 10^, and a blowing rate at the wall of -0.2. Twenty-one nodal points were used 
in the computations. Initial linear profiles for pv(??) were assumed. The accuracy 
criteria were that the computed values of pv(tj), p(??), and a(??) for two consecutive 


(pv)f^ - (pv)^ 1 


(i-lJ 


^ 0.01 for all n 


iterations be within 0.01 at each nodal point (^e.g., 

These accuracy criteria were maintained throughout the entire course of the present study 
of viscous, radiating flows. As indicated by the key in figure 5, the solution converged in 
four iterations which required about 3 seconds of computer time. 


Three constant -density solutions for no blowing are shown in figure 6: an inviscid 
case and two viscous cases for Reynolds numbers of 1 and 10^. The case for Nj^ e = 10^ 
is typical of the Reynolds number of interest, whereas the case for N^ e = 1 is actually 
outside the limits of applicability of the thin-shock-layer equations and is shown merely 
to demonstrate that the present approach yields solutions over the entire range of Reynolds 
numbers. 


The effect of the boundary layer and its extent can be seen in figure 6(a). For 
inviscid flow, there is no boundary layer and a nonzero velocity gradient exists at the 
wall. This gradient, from equation (32), is given by a j = \J(3/Pi< For Nr 6 = 10 5 , most 
of the shock layer is inviscid with only a thin boundary layer present near the wall (to 
rj ~ 0.03), where the velocity gradient slope in y (i.e., da/d??) is a maximum; whereas 
for NR e = 1, the entire shock layer is viscous and no abrupt changes in da/d?? are 
observed. The plot in figure 6(b) indicates that there is little difference in the mass- 
flux distributions for the inviscid case and for the viscous case with Nr 6 = 10^. Also 

-A* 

shown in the figure are the transformed shock-layer thicknesses 6 for the three cases. 
In the inviscid case, the stagnation-line analysis yielded a value for 6 of 0.796, which 
is close to the value correlated by Inouye (6 = 0.78) from inviscid analyses for the entire 
subsonic flow field surrounding a spherical body (ref. 42). For Nr 6 = 10 5 , the trans- 
formed shock-layer thickness 6 is about 3 percent greater than that for the inviscid 
case, whereas for Nr 6 = 1, 6 increased about 33 percent. 

The influence of the blowing rate on the velocity gradient and the shock-layer thick- 
ness y g for a constant density is shown in figure 7 where a is plotted in the physical 
coordinate system. The figure indicates that the inner flow region from the body to the 
stagnation point ^y s - (y)p V= q ~ °» 04 l) is drastically modified as a result of blowing but 
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that the inviscid outer flow is virtually unaffected by the blowing rate. For the constant- 
density model, the shock wave simply moves outward from the body as the blowing rate 
increases, while the character of the inviscid outer flow is independent of the blowing rate. 

The constant -density results demonstrate that for a wide range of blowing rates and 
Reynolds numbers, the solution to the flow equations converges in a minimum number of 
iterations, the nodal-spacing requirements are indeed moderate, and the computational 
times are reasonably short. With this information, the behavior of the overall solution 
with the variable -density iteration can now be examined. 

Nonradiating Air Solutions 

The overall convergence of the flow solution for a variable density is examined in 
this section for equilibrium air without radiation. Solutions for both nonblowing and the 
blowing of air into air are presented. The air model consists of elemental mass frac- 
tions of 0.78 nitrogen and 0.22 oxygen. 

The development of stability criteria is difficult even for the simplest of equations. 
For the coupled system of nonlinear equations which are solved in this study, it becomes 
necessary to rely on the experimental data obtained from the numerical solution in 
order to obtain information about the overall stability and convergence behavior of the 
solution. The results of the information are reflected in the logic presented in the flow 
diagram for the overall solution procedure shown in figure 4. It should be emphasized 
that figure 4 represents the results of the experimental stability study which is used in 
all subsequent studies. The numerical results which led to this particular iterative pro- 
cedure for the variable -density solutions are discussed below. The solutions are for a 
nose radius of 3.43 m and the following free -stream conditions unless otherwise 
noted: 

U’^ = 14.6 km/sec 
p' = 1.6 x 10'^ atm 
p’^ = 2.38 X 10-7 g/ cm 3 

It was observed in the numerical solution that major oscillations occurred in the 
enthalpy and the density, particularly in the viscous region of the flow, as the iterative 
procedure progressed-in-time.— Thes e enthalpy oscillations are shown in figure 8, which 
is a plot of the calculated enthalpies at the first three nodal points" adjacent to the wall - 
and within the shock layer (note that Hj = H w = 0.1 and is constant) versus iteration 
number, where one iteration represents one complete pass through the governing equa- 
tions. From the figure it can be seen that the amplitudes of the oscillations generally 
decrease with i), and although it is not shown here, it was observed that after about five 
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iterations the enthalpy values had converged beyond 77 ~ 0.3. However, the solutions 
near the wall did not converge but merely continued to oscillate even after 30 iterations. 

The density profiles exhibited an oscillatory behavior similar to that of the enthalpy 
profiles but in an inverse fashion (i.e., an overprediction of the enthalpy led 'to an under - 
prediction of the density). This is because there is a strong inverse coupling between the 
density and the enthalpy in the equation of state.' The equation of state for air has been 
correlated by G. Louis Smith (ref. 45) in the form p cc p a b^, where a and b are posi- 
tive exponents. Since p is nearly constant, then an overprediction of the enthalpy 
obtained in the energy equation is a corresponding underprediction of the density from 
the equation of state and vice versa. It is possible to get into a resonant computing mode 
in which the density and the enthalpy oscillate back and forth and can be slowly convergent 
or even divergent. 

From the results it was not clear that the solution was divergent; however, it was 
apparent that if it was converging near the wall, it was prohibitively slow. 

Since the calculations appear to be oscillating about the solution, then by proper 
damping of the calculated quantities it is possible to speed convergence. Fox (ref. 46) 
has discussed the method of underrelaxation, whereby the quantities which are highly 
oscillatory are weighted, or damped, before they are used in subsequent calculations. It 
was found that by weighting the static enthalpy and the density profiles by a certain per- 
centage of their former values and their newly calculated values, rapid convergence of the 
overall solution could be obtained. The weighted values of the static enthalpy and the den- 
sity at each nodal point are computed from the following relations: 

hn = Vi 1 ' 1 * + (1 - 

and 

On = :Vn i ‘ 1) + (1 - d p)Pn’ 

The convergence behavior of the solution for damping factors of 0.5 and 0.9 is 
shown in figure 9. Plotted in this figure are the total enthalpy and density values com- 
puted at nodal point 2. Both solutions converged; however, it is apparent that while over- 
damping will insure convergence of an otherwise oscillatory solution, it can be unduly 
time consuming. For damping factors around 0.2 to 0.3, certain solutions would converge 
but required more iterations than for damping factors of 0.5. For most of the cases 
examined, both with and without blowing, the damping factors of 0.5 appeared to be near 
optimum in terms of generating a converged solution in a minimum number of iterations. 

It was also observed that for the largest blowing rate (-0.2) considered in this analysis, 
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damping factors of 0.7 speeded convergence of the solution by about a factor of 2 over the 
convergence time for damping factors of 0.5. Because of the sensitivity of the computer 
run times to the damping -factors, and since mq attempt was made in the present study to 
optimize the damping factors for a particular solution beyond that which was just dis- 
cussed, subsequent discussions of computer time relate closely to damping factors of 0.5. 
With a thorough study of damping -factor requirements for optimum solutions, run times 
may be improved significantly. 

In order to obtain an indication of the stability of the variable -density solution for 
various blowing rates, three blowing rates, 0, -0.1, and -0.2, were examined. The same 
free-stream conditions given previously and a wall enthalpy value of 0.1 were specified. 
Twenty -one nodal points were used across the shock layer. All solutions converged to 
an absolute density accuracy of 0.01 for damping factors of 0.5. It was observed that the 
number of iterations of the governing equations, and consequently the run times, increased 
with increasing blowing rate. For (pv) w = 0, -0.1, and -0.2, the number of iterations 
required to obtain convergence were 11, 16, and 21, respectively, and the corresponding 
computer times were 1.5, 2.5, and 3.5 minutes, respectively. The computer run times 
are well within reason for computations of the fully viscous shock layer with equilibrium- 
air chemistry. 

The results obtained for (pv) w = 0 and -0.1 are shown in figure 10. On a physical 
scale the boundary layer for (pv) w = 0 occupies about 5 percent of the shock layer, and 
the combined inner layer and boundary layer for (pv) w = -0.1 occupies about 20 percent. 
The results for (pv) w = -0.2 are not shown in the figure. The solution for (pv) w = -0.2 
converged; however, in the inviscid outer region the enthalpy profiles were extremely 
irregular when a central -difference form of the convective term in the energy equation 
was used. Enthalpy values ranging from 0.45 up to about 0.7 (the latter value being 
greater than the free-stream total enthalpy) were observed. To a similar degree, irreg- 
ularities in the density profile and minor inconsistencies in the velocity -gradient profile 
were observed. 

It was observed much later in the study, when the solution to the elemental-diffusion 
equation was required for the study of ablation -products injection, that the energy equa- 
tion and the elemental-diffusion equations are similar and have coefficients on the deriva- 
tive terms which under further dimensional analysis reduce to about the same order of 
magnitude. A^dislTussed _ in _ appendix-B r the elemental-diffusion equation can be unstable 
unless the nodal spacing is made prohibitively small (at) ~ 10“^), and a windward- 
difference formulation is required to obtain meaningful profiles. As a consequence 
of the stability study of the elemental -diffusion equation, a windward-difference form of. . _ 
the energy equation was also used for the radiating cases. 
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Radiating Air Solutions 

The solutions for the viscous, radiating shock layer, including air-to-air injection, 
are discussed in this section. The present results are compared with results from 
existing approaches. 

General results from the present analysis .- The results obtained for (pv) w = 0, 
-0.1, and -0.2 are shown in figure 11 for a nose radius of 3.048 m and the following 
free -stream conditions: 

= 15.25 km/sec 
p'^ = 2.72 X 10’ 7 g/cm3 
p’ = 1.95 x 10 "4 atm 

A wall enthalpy of 0.028, which corresponds to a wall temperature of 3600 K, was 
specified. This temperature is close to the steady-state ablation temperature predicted 
by Smith et al. (ref. 13). Twenty -one nodal points were used for all computations except 
the radiative -flux computations, where 11 points were used to avoid excessive computer 
time. The intermediate values of the radiative heat fluxes required in the solution of 
the energy equation were obtained by linear interpolation. Typical total computer times 
ranged from about 30 minutes for the solution of the problem without blowing to about 
70 minutes for the problem with a blowing rate (pv) w of -0.2. 

A comparison of the radiating solutions shown in figure 11(a) with the nonradiating 
solutions shown in figure 10(a) for slightly different free -stream and wall enthalpies 
shows that the physical shock-displacement distance significantly decreases when radia- 
tion is taken into account. However, the mass -flux distributions were only slightly 
altered in the 77 -coordinate system. The significant decrease in y s when radiation is 
included is due largely to the nonadiabatic effects in the outer inviscid region which 
result in higher values of density (by almost a factor of 2 as shown by comparing fig- 
ures 10(e) and 11(e)) and, to a somewhat lesser extent, the adjustment of the inner flow 
to the wall boundary conditions. On a physical scale the outer inviscid region occupies 
about 95, 75, and 60 percent of the total shock layer for blowing rates of 0, -0.1, and -0.2, 
respectively. 

The radiative heat -flux distributions across the shock layer are shown in fig- 
ure 11(f). Positive heat -flux values indicate net radiative heat transfer away from the 
body; negative values indicate net transfer toward the body. The nondimensional heat 
fluxes at the wall are -0.0434 for (pv) w = 0, -0.0390 for (pv) w = -0.1, and -0.0369 for 
(pv) w = -0.2. These values correspond to dimensional values q^ of -4150, -3740, 
and -3490 watts/cm2, respectively. Thus for air-to-air injection at (pv) w = -0.1 and 
-0.2, the radiative heating rates are reduced 10 and 16 percent, respectively, below the 
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rates for no blowing. The absolute values of heat flux for both blowing rates decrease 
slightly in the inner region in the direction of the body; however, the results indicate that 
the inner air layer is only moderately effective in absorbing the incident radiation from 
the high-temperature outer layer. 

It is observed in figure 11(f) that the radiative heat fluxes reach a minimum near 
the stagnation point, which is expected on the basis of an examination of the governing 
differential equations. Note that if the flow is inviscid, the minimum in the radiative 
heat flux occurs exactly at the stagnation point. This feature provides an important 
self -test of the overall accuracy of the present method, that is, the ability of the method 
to generate thermodynamic properties which, as inputs in the computations of radiative 
heat flux, are sufficiently accurate to yield the minimum in the fluxes near the stagnation 
point. 

More interesting and important checks on the adequacy of the present method for 
the stagnation-line solution can be made by comparison with solutions to the entire sub- 
sonic flow field and by comparison with existing stagnation -line solutions. The first 
exercise serves not only as a check on the accuracy of the solution but also as a check 
of the fundamental assumption that the stagnation-line solution can be decoupled from 
the entire region of influence. 

Comparisons of no -blowing results with existing solutions .- Shown in figure 12 is 
a comparison of the enthalpy profiles generated in the present approach with the 
stagnation-line enthalpy profiles obtained by Suttles (ref. 31) from a solution by a one- 
strip method of integral relations and with those obtained by Falanga and Sullivan (ref. 47) 
from an inverse -method solution. 

The solutions are for no blowing, a nose radius of 3.427 meters, and the following 
free -stream conditions: 

= 14.55 km/sec 
= 2.377 x 10" 7 g/cm3 
p^ = 1.6 x 10‘ 4 atm 

-Both -of the ref erence solut ions are extracted from complete solutions to the entire radiat- 
ing subsonic flow field. The present solution is for viscous-flow , whereas the, re fere nce 
cases are for inviscid flows; however, the comparisons are still meaningful since for 
no blowing the boundary layer -occupies only a small percentage of the shock layer. The 
RATRAP computer code (ref. 10) was used in all three studies. The enthalpy' profile 
from the present solution compares favorably with that from the inverse solution in the 
inviscid region for both N = 11 and N = 21 points (with radiative flux calculations at 
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11 points). The enthalpy profile computed with N = 11 agrees well enough with the 
enthalpy profile computed with N = 21 that the difference is not discernible in figure 12. 

A complete summary of the calculated radiative heat fluxes at the wall is given in 
table 2 for the cases examined in the present study. Also shown in the tabl'e are com- 
parisons with existing sources. Although the inverse solution is considered more accu- 
rate than the integral-relations solution in defining shock-layer profiles, it is interesting 
to note that the radiative heat fluxes at the wall predicted by the three methods fall within 
3 percent of each other. 

Comparisons of air-to-air injection results with existing solutions . - The viscous 
radiating solutions for air-to-air injection with (pv) w = -0.1 are compared in fig- 
ures 13 and 14 for a nose radius of 3.048 meters and the following free -stream conditions: 

= 15.25 km/sec 
pj^ = 2.72 x 10-7 g/ cm 3 
p^ = 1.95 X 10‘ 4 atm 

The wall enthalpy is 0.028, which corresponds to a wall temperature of 3600 K. 

Figure 13 shows results across the entire shock layer for conditions identical with 
those given above. Results from the present analysis are compared with those from 
Rigdon et al. (ref. 12) and Smith et al. (ref. 13) and with unpublished results obtained by 
K. H. Wilson of Lockheed Missiles & Space Company from his method given in refer- 
ence 14. In general, for air-to-air injection, the flow-field results from all approaches 
are in fairly good agreement although there are notable detailed exceptions in the flow- 
field structure which are discussed below. Also, the wall radiative heat fluxes com- 
puted by Wilson, Smith et al., and the present analysis are within 5 percent, whereas the 
results of Rigdon et al. are about 30 percent higher. This disagreement is attributed 
primarily to the differences in the radiation model employed in the first three approaches 
(RATRAP) and that employed by Rigdon et al. (SPECS). 

Although the radiative heat fluxes predicted by Wilson are in good agreement with 
the present method, there was a notable 10-percent difference in the shock-layer thick- 
ness predicted by Wilson and by the other three approaches. The possible sources of 
this difference are associated with the solution of the x-momentum equation and are dis- 
cussed below. 

Tangential velocity gradient: For all the solutions, there is fairly good agreement 
in the tangential velocity gradient in the inviscid outer region and, in general, somewhat 
poorer agreement near the body as shown in figure 13(a). Near the body the velocity 
gradient from the present method is close to that obtained by Rigdon et al. (ref. 12). 
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The major differences near the body in the present solution and in the unpublished results 
of Wilson are due to differences in the computed inner layer and boundary -layer thickness, 
as shown in figure 14(a). 

A possible source of the discrepancy between Wilson's solution for the inner region 

/ a 2\ 

and the present results is the pressure -gradient coefficient [ — £ assumed. The 

\9x2/x=0 

higher the pressure coefficient, the more rapidly the oncoming flow sweeps around the 
body, and thus, the smaller the standoff distance. A coefficient of -3.0 is assumed by 
Wilson, whereas the value of -2.0 is assumed by Rigdon et al. and in the present analysis. 
Smith et al. use a value of -2.5 which has been correlated by Inouye (ref. 42) on the basis 
of inverse flow-field solutions. 

A case was run to examine the effect of the pressure coefficient. The results are 
shown in figure 15. The pressure -gradient coefficient was set to -3.0, and the shock was 
considered to remain concentric. Thus, the velocity gradient behind the shock became, 
by equation (38), a s = N//3/2 = 1.22. The results indicated a decrease in shock-layer 
thickness of 17 percent, which agrees more closely with Wilson's shock-layer thickness. 
However, most of the adjustment occurred in the outer inviscid region (by virtue of the 
higher velocity gradients in this region). As shown in figure 15, there was no significant 
change in the extent of the inner region. 

Two assumptions are made by Wilson (ref. 14) which are not necessary in the pres- 
ent method. First, because of numerical stability problems, Wilson treats the inner 
region as inviscid and begins his fully viscous calculations near the interface between 
the inner region and the boundary layer. A case was run with the assumption that the 
first five nodal points, up to y/y s ~ 0.04, were inviscid. This was accomplished by 
setting p = 0 at these points. There were no discernible differences in the fully vis- 
cous and the inviscid-viscous solutions. This indicates that the assumption of an inviscid 
inner region is valid for the massive -blowing problem. 

The second assumption made by Wilson (ref. 14) in his solution to the governing 
equations is that the product of density and viscosity is constant across the shock layer 
and is given by pp = (pp) w . Since the value pp within the boundary layer for the com- 
parison case examined was about one -half the value at the wall, it was decided to examine 
this assumption"." Two cases~were_run for (pv) w = 0 and -0.1 in which the density- 
viscosity products appearing in the x-momentum equations "and the-energy^equation were 
set equal to -(pp)w at each nodal point. The results of the constant pp cases and the 
corresponding variable pp cases were nearly identical. This is shown in figure 16. 

The results tend to indicate that although viscosity is important in the boundary "layer, -it- - - 
need not be too well defined in the computations since the boundary layer is actually a 
thin transition region which adjusts to the outer and inner flow regions. 
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Again, in figure 13(a), the results from the integral solution of Smith et al. (ref. 13) 
for the total shock-layer thickness agree fairly well with the present results; however, 
their inviscid inner -layer thickness agrees more closely with Wilson's unpublished cal- 
culations, as is seen in figure 14(a). 

Also shown in figure 14(a) are the results obtained by Wilson and Hoshizaki (ref. 28) 
from their integral solution approach. The comparisons indicate that the integral method, 
which was sufficient for radiation studies without blowing, is not adequate for the strong- 
blowing conditions. 

Total enthalpy: The total -enthalpy results across the shock layer are compared in 
figure 13(b). With the exception of the differences in the shock-displacement distances 
which have been noted, the present results compare favorably with the unpublished results 
of Wilson. The enthalpy results of Smith et al. (ref. 13) are fairly good overall. Their 
method is based on a one -strip integral -relations solution for the outer inviscid flow 
coupled with an integral solution for the boundary layer and inner inviscid region. It is 
not expected to yield the detailed flow-field structures obtainable from the other 
approaches (i.e., Wilson (ref. 14), Rigdon et al. (ref. 12), and the present method). 

Temperature: The temperature results obtained from the four approaches agree 
fairly well in the outer inviscid region, but the agreement is somewhat poorer in the 
boundary layer and near the wall, as shown in figure 13(c). Since the radiative heat 
fluxed are strong functions of temperature, the results near the wall are compared in 
more detail in figure 14(b). The present results compare favorably with the unpublished 
results of Wilson. 

The results of the simplified method of Smith et al. (ref. 13) shown in figure 14(b) 
indicate a nearly constant temperature in the inviscid inner region, whereas all other 
solutions indicate an appreciable drop in temperature going toward the body, because of 
radiation exchange. It is unfortunate that because of differences in the radiation models, 
a more direct comparison cannot be made with Rigdon et al. (ref. 12), who also solve the 
fully coupled equations for the viscous, radiating shock layer. However, it is apparent 
from their temperature profile that the radiation model has a strong influence on the 
resulting temperature and predictions of radiative heat fluxes. 

Radiative heat flux: The radiative heat fluxes are shown in figure 13(d). The pres- 
ent results indicate that the inviscid inner air layer is relatively ineffective in reducing 
the radiation to the wall, and as a consequence, the wall heat fluxes of Wilson and of 
Smith et al. (ref. 13), who predict thinner inner layers, agree well with the present 
results. 

A summary of the predicted heat fluxes at the wall for the various air-to-air injec- 
tion rates is given in table 2(a). 
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It is interesting to note that the radiative heat fluxes at the wall converged much 
more rapidly than the detailed thermodynamic properties across the shock layer. For 
example, in figure 17 the heat flux at the wall for (pv) w '= -0.1 -converged wj.thin ±2 per- 
cent of its final value after only five iterations, whereas the density values in the inner 
region (77 = 0.2 and 0.4) were about 50 percent below their final value. This indicates that 
the radiative heat fluxes at the wall for air-to-air injection may not be as sensitive to the 
thermodynamic properties in the inner region as previously expected. If only the gross 
quantity of the heat flux at the wall is of interest, then the density-convergence criterion 
of 1 or 2 percent could be relaxed considerably to provide substantial savings in compu- 
tational time for studies with air-to-air injection. 

As previously noted, the air in the inner region is not too effective in absorbing the 
incident radiation from the high -temperature outer layer. However, since certain chemi- 
cal species of ablators are strong absorbers and emitters, the flux to the wall may be 
sensitive to the relative amounts of these species. If this is the case, then a detailed 
definition of properties across the shock layer is desired. Thus, the density -convergence 
criterion of 2 percent was retained in the subsequent ablation study. 


Solutions of Radiating Flow Field With Ablation Products 

The solutions to the fully coupled, viscous, radiating flow field with injection of 
ablation products are presented in this section. As a result of the stability study (appen- 
dix B) it was decided to retain the windward-difference approximation for this analysis. 
The results of the present analysis are compared with existing solutions. The validity of 
assuming a binary diffusion model is also investigated in this section. 

The heat-shield material is considered to be a carbon -phenolic ablator composed 
of 90 to 95 percent carbon with the remaining elemental mass fractions consisting of 
nitrogen, oxygen, and hydrogen. Exact values of the elemental mass fractions of each of 
the constituents are specified for the particular cases presented. 


Check on assumption of binary diffusion model. - The binary diffusion coefficients 
were computed on the basis of the dominant species present. For the ablator, the domi- 
nant species is atomic carbon, and for air, it is atomic nitrogen; hence, the binary dif- 
fusion-coefficient Jjm diffusion of atomic carbon and atomic nitrogen was used in the 
computations. - 


However, to obtain an indication of the validity of the assumption, two cases were" 
run at identical conditions with the* exception-that diffusion of atomic carbon and atomic 
nitrogen was used in one solution in calculating the diffusion coefficient and diffusion of. _ 
atomic hydrogen and atomic nitrogen was used in the other solution. A comparison of 
ablator mass-fraction profiles obtained from solution of the elemental -diffusion equation 
with the two assumptions is shown in figure 18 for the conditions noted. The two profiles 
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show that the results differ in the mixing region and its extent, with the ablator elements 
extending farther into the shock layer for the hydrogen-nitrogen diffusion coefficient. 

This is to be expected since is almost an order of magnitude larger than D<~:_n, 

as shown in the figure. 

A detailed analysis of the heat -flux distributions for the two cases indicated a maxi- 
mum difference in the radiative heat fluxes of 3 percent, which occurred in the outer 
inviscid mixing region (77 = 0.7), as shown in figure 19. At the wall, the nondimensional 
heat fluxes differed by less than 1 percent (-0.02892 for Djj_n and -0.02917 for Dq_^). 
The shock-layer thicknesses y s were 0.04173 and 0.04168 for Djj_^ and respec- 

tively. Thus, it appears that for the conditions examined in this study, the assumption of 
a binary diffusion model is valid. 

General results from the present analysis .- The results obtained for (pv) w = 0, 
- 0 . 1 , and - 0.2 are shown in figure 20 for a nose radius of 3.048 m and the following free- 
stream conditions: 

U'^ = 15.25 km/sec 
p' rXi = 2.72 X 10 - ^ g/cm^ 
p'^ = 1.95 x 10 - ^ atm 

A wall enthalpy of -0.049, which corresponds to a wall temperature of 3600 K, was 
specified for the cases with injection of carbon -phenolic ablation products with the fol- 
lowing elemental mass fractions: 

(?£ = 0.9207 
= 0.0086 
a Q = 0.0491 
« H = 0.0216 

For the no-blowing case, air is adjacent to the wall, and an enthalpy of 0.028 corresponds 
to the wall temperature of 3600 K. 

In general, the profiles and parameters are much the same as those for air injec- 
tion with the exceptions of the enthalpy profiles for (pv) w = - 0.2 and the radiative heat 
fluxes. 

Even with the windward -difference form of the energy equation, the enthalpy profile 
for a blowing rate (pv) w of -0.2 was irregular near the stagnation point. This solution 
was near convergence when the calculations were terminated on the computer. Appar- 
ently, this blowing rate is near the stability limit for the technique employed herein for 
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the solution to the governing equations. As can be seen in figure 20(f), the density and 
temperature profiles, which are used in the computations of radiative heat flux, are much 
better behaved than the enthalpy profiles. 

It can be seen from the radiative -heat -flux profiles of figure 20(g) that the flux out 
of the shock (at 77 = 1.0) increased with increasing ablator mass injection, whereas for 
air-to-air injection the flux at the shock was nearly insensitive to the blowing rate. (See 
fig. 11(f).) The heat fluxes at the wall, which are of primary interest, decreased with 
blowing rate for the carbon -phenolic ablator; the reductions were 34 and 39 percent of 
the no-blowing radiative heat flux for (pv) w = -0.1 and -0.2, respectively. Thus, the 
ablation products are considerably more effective than air in reducing the heat flux to 
the wall. 


It is not clear whether the large negative heat flux near the stagnation point for 
(pv) w = -0.2 is real or a result of the uncertainty in the calculations of the thermody- 
namic properties. This case is compared with the solution of Rigdon et al. (ref. 12) in 
the following section. 

Comparisons with existing solutions .- Two carbon-phenolic injection cases were 
run for comparison: one at a blowing rate of -0.076 for comparison with the results of 
Smith et al. (ref. 13) and Chin (ref. 2), and one at a blowing rate of -0.2 for comparison 
with the results of Rigdon et al. (ref. 12). The pertinent free -stream and wall conditions 
are noted in figures 21 and 22 for these corresponding cases. In general, the results for 
the thermodynamic and flow properties are in reasonable agreement, but there are notice- 
able exceptions relating to the heat -flux computations which are discussed below. 


It should be noted that the close agreement of the wall heat -flux prediction with 
Chin's computations (ref. 2) may be fortuitous. Chin's earlier radiation model is suffi- 
ciently different from RATRAP to make any quantitative comparison meaningless. This 
is evident from table 2, where it is observed that the predictions of no -blowing radiative 
heat flux from Chin and from the present method differ by 20 percent. 


The comparison of interest for this case is with the results of Smith et al. (ref. 13). 
It is noted in figure 21(e) that the present solution predicts a wall heat flux about 30 per- 
cent lower than that of Smith et al. An important conclusion made in reference 13 is that 
the radiation to the wall is not attenuated as much as previously predicted by Chin and 
-by -Rigdon. etjd._(ref._l %)■ This was attributed to the presence of large percentages of 
CN, which is a strong emitter, in the mixing region! A~comparisonof the-mole fractions - 
of the major chemical species is shown in figure 21(d). The comparisons are good.con- 
sidering the differences in the temperature distributions; however, it is noted that the peak 
^rnole fraction for CN predicted by the present method is only about one -half of that pre- 
dicted by Smith et al. This discrepancy in the level of CN occurs in the region where 
Smith et al. employ a cubic fit to the elemental-diffusion equation in order to join the wall 
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products in the inner shear layer with the air composition beyond the boundary layer. It 
thus appears that contrary to the results of the air-to-air injection study, which did not 
require a very accurate resolution of properties near the wall, the study with ablation- 
products injection indicates that these properties must be well defined near the wall to 
generate accurate predictions of wall radiative heat flux . 

The comparison with Rigdon et al. (ref. 12) indicates reasonably good agreement in 
the radiative heat -flux profiles (fig. 22(a)), although there are significant differences in 
the radiation models employed. The predicted shock-layer thicknesses are within about 
5 percent. A slight waviness is noted in the pv profile of Rigdon et al. (fig. 22(a)) near 
the stagnation point (y/y s = 0.3). In all the cases calculated by the present method, the 
pv profile was smooth in this region. It is not apparent whether this irregularity pre- 
dicted by Rigdon et al. is real or related to the numerical procedure used for- integrating 
out in both directions from the stagnation point. It is noted that their solution'also exhibits 
an irregularity in the temperature profile near the stagnation point as does the present 
solution. 


Summary of Wall Radiative Heat Fluxes 

A summary of the wall radiative heat fluxes calculated by the present method for 
the various conditions examined in this study and corresponding predictions from previous 
investigations are given in table 2 for comparison. With the exception of the no-blowing 
results of Suttles (ref. 31) and the air-injection case of Wilson (unpublished), the present 
method predicts heating rates slightly to considerably lower than those calculated by 
previous investigators. 

These results are also summarized in figure 23, which illustrates the effectiveness 
of the carbon -phenolic ablator products in reducing the wall radiative heat fluxes. As can 
be seen from the figure, there is general agreement of all sources that air injection is 
moderately effective in reducing the heat transfer to the wall and, with the exception of 
the results of Smith et al. (ref. 13), that the ablation products of the carbon-phenolic heat 
shield are highly effective in reducing the fluxes. 

As can also be seen from figure 23 and table 2, there is no consistent agreement in 
any of the sources when both the no-blowing and blowing results are considered. For 
sources which use different radiation -transport models (Rigdon et al. (ref. 12) and Chin 
(ref. 2)), this disagreement is due, at least in part, to the differences in the radiation 
models. For sources which use the same model (Smith et al. and the present method), 
this disagreement is attributed to differences in the temperature distributions near the 
wall. The data available were insufficient to evaluate fully the unpublished results of 
Wilson, with the exception of the assessment of the validity of the assumptions required 
in his analysis. These assumptions have been previously discussed and shown to be valid 
for air-to-air injection with (pv) w = -0.1. 
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For air-to-air injection, the present results indicate a 10-percent reduction in the 
no-blowing wall radiative heat flux for (pv) w = -0.1 and a 16 -percent reduction for 
(pv) w = -0.2. For injection of carbon-phenolic ablation products at corresponding blow- 
ing rates, the reduction is 34 and 39 percent, respectively. Since ablation rates are 
typically expected to be about 0.1 (see Smith et al. (ref. 13)), then the 34-percent reduc- 
tion in the radiative heat flux represents a significant savings in heat-shield weight. 

CONCLUDING REMARKS 

An implicit finite -difference scheme is developed for the solution of the fully coupled 
equations for the viscous, radiating stagnation streamline with the effects of strong blowing 
included. Solutions are presented for both air-to-air injection and ablation -products injec- 
tion at rates up to 20 percent of the free-stream mass-flow rate. The free-stream condi- 
tions examined are typical of interplanetary return conditions into earth's atmosphere near 
the point of peak radiative heating in the entry trajectory. A detailed radiative -transport 
computer code (RATRAP) which accounts for both continuum and line radiation exchange 
processes is utilized in the study. 

With a minimum number of assumptions for the initially unknown parameters and 
profile distributions, convergent solutions to the full stagnation-line equations are rapidly 
obtained by a method of successive approximations. No singularities exist in this formu- 
lation of the finite -difference equations. Damping of selected profiles is required to aid 
convergence of the massive -blowing cases; however, even for these cases, no "patching" 
of the viscous and inviscid regions is required. The results demonstrate that windward 
differencing of the convective term in the elemental-diffusion equation is required for a 
stable solution to this equation. Although the central-difference approximation to the 
energy equation yields satisfactory solutions when radiation is included, the results are 
considerably improved for a blowing rate of 20 percent of the free-stream rate when the 
convective term is windward differenced in this equation also. 

Comparisons are made with currently existing solutions to the equations for the 
radiating shock layer. The present method predicts lower wall radiative heat fluxes for 
carbon -phenolic ablation than those predicted by previous investigators. 

The-resultS-indicate that the ablation products are highly effective in blocking the 
incident radiation from the high -temperature "outer - layer of the-shock.^ Eor_blowing rates 
of 0.1 and 0.2, typical reductions in radiative heat flux at the wall range from 34 to 39 per- 
o cent of the values for no blowing. - . _ 

For air injection, the inner air layer is shown to be relatively ineffective in blocking 
the incident radiation; hence, the thermodynamic properties need not be as well defined for 
air injection as for ablation -products injection. 
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The validity of assuming a binary diffusion model was examined in the present 
analysis, and although a study with a multicomponent diffusion model may remain mean- 
ingful for conditions and ablators other than those examined in the present study, the 
results indicate that the multicomponent diffusion model is not required. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., July 28, 1972. 
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APPENDIX A 


RADIATION MODEL 

The RATRAP radiative -transport model developed by Wilson (ref. 10) has been 
discussed in great detail by Wilson (ref. 10), Suttles (refs. 37 and 31), and Wilson and 
Hoshizaki (ref. 28). A summary of Suttles' discussion on the development of the govern- 
ing radiative -transport equation employed by RATRAP and some details of the radiation 
model are given in this appendix. 

The RATRAP code, which includes the detailed line radiative calculations (some- 
times referred to as RATRAP II, since it has been upgraded by its originators, see Wilson 
and Hoshizaki (ref. 28)), is written for the calculation of the radiative heat flux at any point 
within a planar (tangent) slab in which the thermodynamic properties vary in only the 
direction normal to the slab. Local thermodynamic and chemical equilibrium is assumed. 
A distribution of two thermodynamic variables plus the elemental mass fractions for gas- 
eous mixtures of carbon, nitrogen, oxygen, and hydrogen is required for the radiation 
computations. 

In order to evaluate equation (17) for the radiative heat flux, it is necessary to cal- 
culate the spectral linear absorption coefficients a v , which are functions of the thermo- 
dynamic properties (p and T) and the number densities of the chemical species. Twenty 
chemical species are considered in the thermodynamic calculations. They are C2, N2, 

0 2 , H 2 , C, N, O, H, CO, CN, C 2 H, C3H, C4H, HCN, C 2 H 2 , C~, C + , N + , 0 + , and H + . 

The total spectral absorption coefficient is separated into continuum and line con- 
tributions in the RATRAP calculations. The continuum spectral absorption processes 
considered are the free-free transitions (acceleration of free electrons in the vicinity of 
atoms and ions) of C, C + , N, N + , O, 0 + and the free-bound and bound-free transitions 
(deionization and ionization, respectively, of the H atoms and positively ionized particles). 
Molecular bands are also included in the continuum calculations. These molecule con- 
tributions and frequency intervals from Wilson and Hoshizaki (ref. 28) are tabulated 
below. 


Molecule 

contribution 

Photon energy range Jtfv for significant 
absorption, eV 

— Hg-Werner 

11 15.494 

Photoionization 

15;494-s.#„_<_25_ 

C 2 Swan 

1.8 ^ ^ 6.0 

Fox-Herzberg 

1.8 S = 5.35 

Mulliken- - - 

5.35 S $ v j'6.0 

Freymark 

1.8 2 Hi/ s 6;0 - - _ 

CN violet 

2.0 2tfv 2 6.0 

CO 4th positive 

O 

VII 

*6. 

VII 

c— 

N 2 Birge-Hopfield 

11 S s 14.2 

0 2 Schumann -Runge 

7 £ K* = 9.2 
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APPENDIX A - Concluded 


For the calculation of the continuum contribution to the heat flux, the absorption 
coefficients for the individual species are weighted with their respective number densities 
and approximated by curve fits over frequency ranges. The heat -flux equation (17) is 
subsequently evaluated by numerically integrating over y and over the frequency, with 
11 values of y and 31 values of frequency used. 

The atomic -line radiation component arises as a result of the many bound electronic 
transitions which occur in the atomic nitrogen and oxygen species (line radiation from 
carbon and hydrogen species is not included in RATRAP). This component is obtained by 
first grouping certain line contributions within various frequency intervals. The net line 
radiation is then calculated by summing the contributions from these line groups. Eigh- 
teen line groups, with a total of 65 lines, are used in RATRAP. 

The total radiative heat flux at each point is obtained by adding the continuum and 
line radiation contributions. 
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APPENDIX B 


STABILITY STUDY OF THE ELEMENTAL-DIFFUSION 
AND ENERGY EQUATIONS 


The solution to the elemental -diffusion equation is required in the analysis of the 
shock layer with ablation-products injection. Use of the central -difference form of the 
diffusion equation results in a numerical instability. This numerical instability is a 
strong oscillation across the shock layer in the computed elemental profiles. This insta- 
bility is eliminated by windward differencing of the convective term in the diffusion equa- 
tion. This appendix presents a stability analysis of the two formulations of the diffusion 
equation and also the energy equation, which is analogous to the diffusion equation, with 
an additional term for radiation. 

The properties vary rapidly across the shock layer. In order to obtain some insight 
into the stability of the diffusion equation and the associated nodal -spacing requirements, 
the. elemental -diffusion equation (35) is simplified by assuming k^p^Dj^ and 6/c^pv to 
be constant, so that equation (35) simplifies to 

d% + Jv_cg =0 (B1) 

d7j2 pDj2 d?7 

Equation (Bl) is cast in finite -difference form by using central differences: 


a . - 2 a + oi i ? ol i - <x 1 
n+1 n n-1 5v n+1 n-1 _ q 


M 


pD 12 2 At? 


(B2) 


For - v A? L = Constant, this equation may be solved by assuming the solution 
2pDf2 

. n 


“n = CX 


(B3) 


Substitution of equation (B3) into equation (B2) and simplification results in the charac- 
teristic equation for X: 


i ~ + hA j\x 2 - 2x + 1 1 - hAl) = 


2pDi2/ 
The- solutions for X are 


2pDi. 2/ 


0 


(B4) 


~\ 




x 2 = 


1 _ 6v A?7 
2pDj 2 

1 , 5v At? 
2pD i2 ^ 




(B5) 
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APPENDIX B - Continued 


The Xj solution corresponds to a constant error across the shock layer and, as such, 
is removed by the boundary conditions. The X 2 solution is an oscillation if X 2 < 0, 


which implies that 


5v At? 
2pDi2 


> 1. 


If X 2 > 0, this contribution has a monotonic variation 

across the shock layer and is thus eliminated by the boundary conditions. The condition 
for eliminating the undesirable oscillation is thus 


pv5 A 77 
2p 2 Di2 


< 1 


(B6) 


and 


The nodal-spacing requirement for stability can be established by noting that 
6 = 0[l] 

pv = 0[l] 


In this paper, a typical value of p is 20, and a typical value of D^2 is 10"®, so that 

2p^Df2 1 - _3l 

A77 < ^ = O 0.8 x 10 6 

pv 6 L J 


This nodal spacing corresponds to approximately 1250 nodal points. 

A case was run in which 1001 nodal points were used for the solution of the 
elemental -diffusion equation, with the same 21 points used for the other equations. 
Although some improvement was noted in the solution, the instability was not eliminated. 

The numerical results which substantiate this stability analysis are shown in fig- 

-R 

ure 24. The computed mass fractions for a diffusion coefficient of 10 are shown in 
figure 24(a) for a nodal spacing of 0.05. The solution is clearly unstable. The results 
obtained when the coefficient was increased to 1.25 x 10'^ and At] was set to 0.1 (which 
accomplishes the same stability effect as decreasing the nodal spacing to 4 x 10 _ 4) is 
shown in figure 24(b). The sharp oscillations have vanished, but aTp still exceeds the 
value of 1.0. The case was rerun for a nodal spacing of 0.02 with no noticeable improve- 
ment. Only when Dj2 was raised to 10 -2 did c?p, monotonically decrease from 1.0 at 
the wall to 0 at the shock. 


Since the central -difference form of the elemental -diffusion equation was inadequate, 
a windward -difference form for the convective term was employed. It has been shown 
that windward differencing of the convective term stabilizes certain fluid-flow computa- 
tions. (See, e.g., ref. 48 or ref. 49.) 

Equation (Bl), when written in finite -difference form with windward differencing of 
the convective term, becomes 
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APPENDIX B - Continued 


a i - 2 o' + a i 
n+l n n -1 


“n+l - “n 


(A 17) 2 


Sv 

P D 12 Ar? 


= 0 


(B7) 


In equation (B7), it is assumed that v is positive, so that a fluid particle would flow 
from the (n+l)th node to the nth node. If v is negative, the appropriate windward- 
difference form is 


a- , 1 - 2a + a < 
n+l n n -1 


6 v 


a 


n 


a 


n -1 


(At ?) 2 


P d 12 At? 


= 0 


(B 8 ) 


As before, a solution form (B3) is assumed, and results in the characteristic equation 
for X: 


+ £v6_At Z \ a 2 _ / 2 ± £v6_^V 1 = 0 
P 2 ° 12/ \ p2D 12/ 


(B9) 


The plus sign corresponds to positive v when equation (B7) is used, and the minus sign 
corresponds to negative v when equation (B 8 ) is used. The solutions for X are 


Xi = ! 


x 2 - 


1 ± 


pv 5 At? 
p 2 e>i 2 


(BIO) 


Because both roots of X are positive, there are no oscillations in the solution when 
windward differencing of the convective term is employed. 

A physical explanation given for using windward differencing of the convective term 
is that if fluid is flowing downstream, the upstream cell influences the downstream cell 
more than the downstream cell influences the upstream cell. The mathematical reason 
for windward differencing is that higher order harmonics (ref. 48) that are introduced into 
the solution because of finite -difference approximations to the differential equations decay 
exponentially, whereas for central -difference formulations these higher order distur- 
bances can be exponentially amplified. 

The'results- obtained for both^ central- and windward-difference forms of the diffu- 
sion equation are shown in figure 25. The computations wereTor a viscous, radiating 
shock layer, and a binary diffusion coefficient for atomic hydrogen and atomic nitrogen 
was assumed. The windward -difference formulation yields a stable solution in which 
Oj, monotonically decreases across the shock layer. " ------- 

A check case was run to establish the effect on the o’-, profiles of the damping 
that is introduced by the windward-difference formulation. This damping, which is a 
form of artifical viscosity, tends to smooth, or smear out, the gradients more than is 
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APPENDIX B - Continued 


natural. Both central- and windward -difference solutions were obtained for diffusion 
coefficients which yield stable solutions for central-difference approximations. The 
results from these two formulations, shown in figure 26, indicate that the windward- 
difference formulation does not introduce any noticeable artificial viscosity for this case. 


When constant density and constant transport properties are assumed, the energy 
and the elemental-diffusion equations (eqs. (34) and (35)), respectively, can be written as 
(where the convective term on the right-hand side of the energy equation has been 
neglected) 

p .,dH_ PM d 2 H 
dT? ® N Re N Pr dr?2 


and 


da-p 

pv — - 
drj 


p 2 D 12 d 2 Q! F 
6 dry 2 


PM 

S N Re N Sc 


d 2 <* F 

drjZ 


where Ng c , the Schmidt number, is given by 


jj = M 
iN Sc 


p’D 


12 


Since the Prandtl and Schmidt numbers are typically of order unity, then the energy equa- 
tion and the momentum equation are controlled by the same stability requirements as the 
elemental-diffusion equation. When radiation becomes a significant contributor to energy 
transport, it acts as an additional mode of conduction. Thus, in central-difference form, 
radiation increases the stability of the energy equation beyond that of the elemental - 
diffusion equation, so that a well-behaved diffusion solution assures a well-behaved energy 
solution. 


As a result of this analysis, the convective term in the energy equation f pv ^2) was 

V d? 7/ 

recast in windward -difference form and the nonradiating equilibrium -air solution for 
(pv) w = -0*2 was rerun. The profiles of enthalpy, density, and velocity gradient became 
smooth. The results of the solution for the x-momentum and energy equations are plotted 
in figure 27 and compared with the previous central -difference solutions. 


A significant change in the shock layer thickness also occurred when the windward - 
difference form of the energy equation was used to improve the profiles. The shock- 
layer thickness y s was reduced by about 12 percent from 0.0703 to 0.0622 in going from 
central- to windward -difference form. It was found that this reduction in y s was due 
primarily to the large blowing rate which introduced the strong oscillations in the profiles 
for the central -difference formulation. The nonradiating cases for (pv) w = 0 and 
(pv) w = -0.1 were rerun with the improved formulation of the energy equation. There 
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APPENDIX B - Concluded 


was no change in y g for (pv) w = 0, and there was a 3 -percent reduction in y s for 
(pv)w = -0.1 ...... 

The strong oscillations in the enthalpy profiles were not present in the radiating 

solutions. However, the (pv) w = -0.1 case for = 15.25 km/sec was rerun to 
determine the effect of the windward-difference approximation in energy equations. The 
radiative flux at the wall was unchanged and the shock standoff distance increased by only 
2 percent. 

As a result of this stability study it was decided to use the windward -difference 
approximations to the convective terms in the energy and elemental-diffusion equations 
for the analysis with radiation and injection of ablation products. 
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TABLE 1.- MOLECULAR CONSTANTS FOR 
DIFFUSION COEFFICIENT 


Species 

Molecular constants 

£ !/k* ; K 

f A 

a v A 

Mi, g/ g mole 

C 

30.6 

3.385 

12.01 

N 

71.4 

3.298 

14.01 

H 

37.0 

2.708 

1.008 




TABLE 2.- SUMMARY OF WALL RADIATIVE HEAT FLUXES 

(a) Air injection 


(pv) w 

Source 

q R,w’ watts /cm 2 

u'oo = 14.55 km/sec; p'^ = 2.377 X 10 " 7 g/cm 3 ; = 1.6 X 10' 4 atm; R^ = 3.427 m 

0 

Suttles (ref. 31) 

-2860 


Falanga and Sullivan (ref. 47) 

-3000 


Present analysis 

1 -2950 

UU = 15.25 km/sec; p’^ = 2.72 X 10“ 7 g/cm 3 ; = 1.95 X 10‘ 4 atm; = 3.048 m 

0 

Rigdon et al. (ref. 12) 

-5990 


Smith et al. (ref. 13) 

-4170 


Present analysis 

-4150 

-0.1 

Rigdon et al. (ref. 12) 

-5130 


Smith et al. (ref. 13) 

-3991 


Wilson (unpublished) 

-3660 


Present analysis 

-3740 

-0.2 

Rigdon et al. (ref. 12) 

-4630 


Present analysis 

-3490 


(b) Carbon-phenolic injection 


(pv) w 

Source 

q R,w’ watts /cm 2 

U’oo = 15.25 km/sec; 

Poo = 1-77 X 10 -7 g/cm 3 ; p’^ = 

1.21 x 10- 4 atm; Rj^ = 2.56 m 

0 

Chin (ref. 2) 

-3030 


Smith et al. (ref. 13) 

-2795 


Present analysis 

-2540 

-0.076 

Chin (ref. 2) 

-1642 


Smith et al. (ref. 13) 

-2188 


Present analysis 

-1620 

Uoo — 15.25 km/sec; | 

- 2 - 72 * iO' 7 g/cm 3 ; p^ = 

1.95 X 10" 4 atm; R^ = 3.048 m 

0 

Rigdon et al. (ref. 12) 

-5990 


Smith et al. (ref. 13) 

-4170 


Present analysis 

-4150 

-0.1 

Present analysis 

-2790 

-0.2 

Rigdon et al. (ref. 12) 

-2670 


Smith et al. (ref. 13) 

-3246 


Present analysis 

-2560 
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Figure 2.- Flow -field coordinate system in the transformed coordinates. 



Figure 3.- Finite -difference network. (Nodal points are equally spaced in 77 .) 
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Figure 4.- Flow diagram of the overall solution procedure 
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(a) Velocity -gradient distributions. 



(b) Mass -flux distributions. 

Figure 6.- Constant-density solutions for inviscid and viscous flows 

without blowing, p = 20 * 
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Fjgure 9.- Effect of profile damping on the convergence of the variable -density solution. 

Equilibrium air without radiation. 
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(b) y -momentum equation. 
Figure 10.- Continued. 
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Figure 11.- Equilibrium -air solutions with radiation included. 




n (P v ) w * 0; 6 = 0.858; y s ■ 0.0334 

O (pv) w » -0.1; 6 ■ 1.629; y s = 0.0408 

A (pv) w » -0.2; 5 m 2.544; y g » 0.0456 

Solid symbols indicate approximate 
location of the stagnation point 
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(d) Energy equation. (Central -difference form of convective term in energy 
equation used to generate these solutions.) 

Figure 11.- Continued. 




(e) Equati 
Figure 11. 


(pv) w ■ 0; 5 = 0.858; y s = 0.0334 

(pv) w = -0.1; 6 = 1.629; y s * 0.0408 

(pv) w = -0.2; 6 = 2.544; y s > 0.0456 

Solid symbols indicate approximate 
location of the stagnation point 
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(a) Tangential velocity gradient. 

Figure 13.- Comparison of equilibrium -air solutions for radiation and 

air-to-air injection. 
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Figure 13.- Continued. 
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Figure 17.- 


(a) Density behavior. 
Study of convergence of density 


and radiative heat flux 
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— Present method y g = 0.0404 

Smith et al. (ref. 13) y g = 0.0431 

O Chin (ref. 2) y • 0.0410 
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(e) Radiative heat flux. 
Figure 21.- Concluded. 
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(b) Energy equation. 

Figure 27.- Concluded. 
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